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I INTRODUCTION 


Digital signal processing is a field which 1s rapidly expanding due to advances in 
modern technology. Essential to this field are digital filters. Modeling these filters con- 
stitutes much of the effort involved in digital signal processing. The filters provide a 
transfer function which describes the relationship between filter input and output. 
Autoregressive (AR), moving average (MA) and combination autoregressive moving 
average (ARMA) models are widely used to represent the transfer function of a digital 
filter. A filter transfer function is commonlv described in direct form. This form is a ratio 
of two polyvnomuals, usually of the form, 

B(z) by thy zo tb, 


A) =——- = OO (ie) 
A(z) Weve ae 





The above equation describes an ARMA model of order (s,t) where s is the order of the 
denonunator and t the order of the numerator. The a, parameters form the 
autoregressive portion of the ARMA model. The 6, parameters form the moving average 
portion of the ARMA model. [fall the autoregressive parameters are zero, then the filter 
transfer function H(z) is strictly a moving average process of order t. If all the moving 
average parameters are zero except for >) equal to one, then the filter transfer function 


1s Strictly an autoregressive process of order s. 


A. OBJECTIVES OF THE THESIS 

Fundamental to the design of digital filters is estimation of AR, MA or ARMA 
parameters. Accurate and efficient parameter estimation has been the subject in much 
of the related digital signal processing literature [Refs. 1 ,2.3]. The first objective of this 
thesis 1s to confirm the proposed ARMA parameter estimation algorithm of [Ref. 4: pp. 
619-621], which leads to the design of a new ARMA digital lattice filter. The proposed 
algorithm 1s a generalization of the Mullis-Roberts criterion for parameter estimation 
known as the modified least squares problem [Ref. 5: pp. 227-228]. The algorithm uses 
two recursive formula to estimate the parameters. One is an AR recursive formula which 
estimates ARMA parametes as the AR order is increased bv one. The other is an MA 
recursive formula which estimates ARMA parameters as the MA order is increased by 


one. [This algorithm is unique in that it allows for the design of an ARMA model with 


arbitrary AR and MA orders with no dependency of an AR model on an MA model or 
vice versa. The ARMA digital filter designed from the proposed ARMA parameter es- 
timation algorithm is in the form of a lattice structure. Lattice realizations of filters are 
widely used and provide excellent analysis of prediction errors [Refs. 6: pp. 165-168.7]. 
Grav and Markel developed an algorithm which produces lattice realizations of filters 
{rom the direct 1ormn | ice. 0) 

The second objective of the thesis is to make the proposed ARMA digital lattice 
filter of [Ref. 4: p. 662] adaptive. An adaptive lattice filter is one in which the lattice 
coefficients are automatically adjusted by an adaptive algorithm to yield the optimum 
filter design. The adaptive lattice algorithm derived in this thesis is based on the widely 
used least mean square (LMS) algorithm. Adaptive filters have many applications [Ref 
9: pp. 7-31] including. 

[2] Sssteni identification: 


2. Digital representation of speech. 


oo) 


Adaptive auotoregressive spectrum analvsis. 


4. 


Adaptive detection of a signal in noise of unknown statistics. 


Eeho cancellation. 


now 


Adaptive line enhancment. 


~~] 


Adaptive beamforming. 


The need for an adaptive filter is made apparent bv considering a filter of fixed design 
Which is optnuzed for given input conditions. In practice, the complete range of input 
conditions may not be known or could change from time to time. A filter of fixed design 
would not produce optimum results under these conditions. An adaptive filter. which 
yields optimum results given changing input conditions, will give superior performance 
to one of fixed design. 

The last objective is to analyze convergence properties of the derived adaptive lattice 
algorithm. This is accomplished by computer implementation of the adaptive algorithm. 
The output of a known transfer function is compared to the output of the adaptive 
lattice filter given a common input. Plots of the error between the two outputs and 


lattice coefficient convergence are obtained. 


B. ORGANIZATION OF THE THESIS 
Chapter II is designed to present the ARMA parameter estimation algorithm and 


ARMA digital lattice filter proposed in [Ref. 4: pp. 617-628]. Computer simulation of the 


tJ 


algorithm was performed and results are shown. A brief review of the Mullis-Roberts 
criterion is provided to establish a reference for expanding this criterion to the proposed 
ARMA parameter estimation algorithm. An adaptive lattice algorithm is derived in 
Chapter II] which makes the proposed ARMA digital lattice filter adaptive. The adap- 
tive lattice algorithm is efficient and accurate. Chapter IV contains experimental results 
which show convergence aspects of the adaptive lattice algorithm when applied to 
ARMA lattice filters. Conclusions about the proposed ARMA parameter estimation 


algorithm as Well as the derived adaptive algorithm are discussed. 


Il. ARMA DIGITAL LATTICE FILTER 


In this chapter, we will review the Mullis-Roberts criterion for solving linear ap- 
proximation problems and introduce analysis equations of the ARMA digital lattice fil- 
ter. The criterion used in the formulation of the ARMA digital lattice filter is a 
generalized form of the Mullis-Roberts criterion [Ref. 5: pp. 227-228], which has been 


given as a modified least mean square problem for ARMA parameter estimation. 


A. MULLIS-ROBERTS CRITERION 

The Mullis-Roberts criterion evolved from considering second order statistics in 
conjuction with first order information about a process to obtain filter approximations. 
Consider the bounded impulse response sequence A = {/p, h,, ... } containing [irstyomeas 


information about the filter # having a frequency response function, 


He N= Dae. (2.1) 
k=0 


Let { u, } be a zero-mean, unit-variance, white-noise sequence and { y, } be the output 
process corresponding to the input { uw, }. then we have the following convolutional re- 


lationship eiven by. 


Polen (292) 
0 


Vie 


Ud 


la 


k 


Second order information about the filter /& 1s obtained from the autocorrelation se- 


quence { r, } given as 


se Eye) ee ae (2aey 
=) 


From equations (2.1) and (2.2), the second order interpolation problem 1s to find a 
lowest order recursive filter which matches the data { fy, ... , /i,. %,--. «%, } ame lem amma 


presents the first order information and r, the second order statistics. 


Let us now consider the case where only first order information about a process 1s 
Palen ae nat 1s, Ciel) an timpuise response sequence { /,, 1,.... } , we Want to find a re- 


cursive filter of the form, 


7 


=I = 
Go FZ IE AZ : 


vn » —k 
f(z) i hyz =} = 
oe hs sn wea male par 


k=0 


n 


which approximates H(z) and therefore the impulse response sequence { /,. fy, ... } . We 
also desire the frequency response H(e“) to approximate the desired response //(e”). 


Suppose that H(e”) 1s chosen such that it minimizes the integral squared error, 


el) — He)Pdo = hf? 2. 


tJ 
Cay 
wee” 


) 


de fh 


Using the Parseval relation, we can obtain an alternative definition of the approximation 


error, 


oO 


— | Oe ee cs = one = ) (a2, (2.6) 


k=O 
Ii the filters H(z) and H(z) are driven by the same white noise source, equation (2.6) 
describes the output error between the filters which we write as, 


Wa-Al? = £G,-fY = Ee) (2.7) 


where y, and }, are the outputs of the respective filters when driven by the same white 
noise source as in (2.2). Minimizing (2.7) is a nonlinear programming problem requiring 
the entire impulse response sequence. As a result, computational efforts for obtaining a 
solution are inefficient. A modification to the problem was introduced [Ref. 5: pp. 
227-228] which considered a cost function that is quadratic in the coefficients of the re- 


cursive filter given by (2.4). The modification is described as follows. Let 


N =" - 
Oe Gag ce) (2.8) 


ane 


Az) = =‘ + ao +e+az") (29) 


be the numerator and denominator polynomials, respectively. in (2.4), where .\ = max 


(n1.7). The task now is to find coefficients which mininuze the quadratic form. 


Ee) =s— | [ey Ale) — Ole Pde (2.10) 


This 1s a standard quadratic mininuzation problem whose integral can be expressed in 
terms of the coefficients of polynomials A(z) and Q(z) and the yi@ams 


{Ho es Ags Now oe Yn be TElating to the filter //(z). This problem 1s shown i Figiresiaaas 


equation (2.10) 1s known as the Mullis-Roberts criterion. 





Figure 1.  Mfodified least squares problem 


B. GENERALIZED MULLIS-ROBERTS CRITERION 


In order to define the new criterion used for ARMA parameter estimation we con- 


sider the following transfer function with input sequence { x(k), k = 1,2,... } and output 
SEQUENCE -yiA). A 0102 a eer cice 
H( a (2.11) 
Th) 





Figure 2. Equivalent input/output model 


where /7,(z) and H,(z) are reference polvnomuals which we desire to model. An equivalent 
model 1s shown 1n Figure 2 where uw(A) 1s an intermediate signal, and the realization 1s 
Similar to that of the direct form realization II [Ref. 10: p. 151]. Let the intermediate 
sequence { wk), k =1,2.... } be a zero-mean white gaussian process. The model of Fig- 
ure 2 can then be transformed into the model of Figure 3 with u(k) as the common input 


to both transfer functions //,(z) and //,(z). 





Figure 3. Transformed model 


Note that we earher defined transfer functions //,(z) and H,(z) as polvnomuals of the 


reference model. For each transfer function an estimation polynomial is defined such 


That 


A(z) = 1+ Crea tetas estimates H(z) (2. 
B(z) = by + byz7 to +52! estimates He) (2.1 
where a, and 6, (¢=1.2,...,s) and Y=1,2.....2) , are the AR and NIA parameter eee 


specuvely. of the combined ARMA model formed by A(z) and B(z). This refined least 
squares problem 1s shown in Figure 4 and is a generalized form of the Mullis-Roberts 
criterion. If the reference model polynomial //,(z) in Figure 4 is equal to unity. we obtain 
the Mullis-Roberts criterion shown in Figure 1. Therefore. bv including releremes 


polvnomual //,(z), the new criterion for ARMA parameter estimation becomes, 


2 4 
Cuma a ate 
Ex. = Eo FIHOAG)— HBS 2.14) 





Ficure 4. Refined least squares problem 


Minimizing E,, 1s accomplished by calculating the coefficients of A(z) and B(z) which 
minimize e(k)? of Figure 4. Another form of eq (2.14) 1s obtained bv applying Parseval’s 


theorem and 1s expressed as, 
E,. = EL(A(z)s(&) — Blz)x())] (2.15) 


Which 1s obvious from Figure 4. The coeflicients a, (j= 139.) mance (j= 1, 


Which nuninuze £,, can then be calculated using the normal equations for ARMA 


5,f 


parameter estimation. In order to obtain the normal equations for the problem in (2.15). 


let us define the following: 


a,,=[a,..a,] and b,,=[d, ... 4) (26) 
are the vectors of AR and MA parameters, respectively, 
er Cte eK S) eK) 23s c(h —'t) 4 tT) 
is the data vector consisting of both input and output data elements and 
R, ,= Elh, ,(4)’ hy (4)] (2.18) 
is the data autocorrelation matrix. The criterion 1s to minimize the mean squared error 


Eg = EL c(h) J= EL LL 1 ase & Doe Jee P| 


_ (2.19) 
= el | (A) a [ as bo bs: ] Hey ig 


Now following the standard calculus of variation optimization procedure for minimizing 


E.., [Ref. 11: pp. 100-110], vields the normal equations in the matrix form, 


See oe, | Ro — menue, 0-00} (2.20) 


It is interesting to note that if E,, in equation (2.14) is zero, then we have the following 


equality, 
H(z) =~ = —— OLD 


Sertie estimate for the total reference mode] H(z) is the ratio of B(z), the MA part and 
A(z), the AR part of an estrmated ARMA model. 


C. ARMA PARAMETER ESTIMATION 

Now that the criterion for ARMA parameter estimation has been established, the 
solution method to estimate the ARMA model parameters to minimize equation (2.14) 
is considered. Let x(k) and }(k) be the input and output signals, respectively, of the es- 


timated ARMA model. Using a difference equation representation, this process is de- 


scribed by 


s 


y(k) = ~) 4 (A +) x x(k — 1) (2.22) 


= |] ra 


~, 


For these input and output signals we define four estimation, or prediction, models as 


follows. The forward estimation signal for x(A) is defined as, 


-_——s 
ty 
tv 
(oP) 

ae” 


t 5 
xk) =— ») b* x(k — i) + »y a; VA 
| 


i=] 


where 67 and a’ are the corresponding estimation parameters. The forward esitmation 


signal for (xk) 1s similarly defined as, 


Ay f 
==) af Hk-D+ DoF xh=2 (2.24) 
i=] 


f= 


The backward estimation errors for x(k — 1) and y(k — s) are then given by, 


~] —~] 

en) ae (ei 5 a? (k —J) (2.25) 
i=0 j/=0 
s—] —1 

(A — S)=— » ay JA -yJ + y b? x (K — 1) (2.26) 
J=0 =o 


Where the superscripts g and d indicate the backward estimation parameters for x and 
Male DeCucin 

From these estimation models, we can now obtain the four prediction errors at any 
given ume k. These errors are expressed as differences between the predicted value and 


actual value of an input or output signal, namely. 


ye (k) = —x(k) + Xk) (2.27) 
vs k) = 3k) — yA) (2.28) 
8.,(k) = x(k — 1) — (k= 1) (2.29) 
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d, (k) = jk — 5) — Yolk — 5) (2.30) 


We now use the vector notation to simplify the expressions for prediction errors. In the 
following. the forward error elements corresponding to both x and y form a vector 


v,(k), given by 
V.AN=L —v tk) 2K =h, (dC, (2.31) 
and the backward error vectors are given by 
g, {k) =—h, (A) Gi, (2.32) 
d, (k) =h, (k) Df, (2.33) 


where C,, and G,, and D,, are the forward estimation parameter matrix and backward 


estimation parameter vectors, respectively, defined as 


z (Mae fa el be... by (2.34) 
eae ic | Ombe 2c . 
Go SLC es eo (225) 
De ice ae, £ to. 62, 0] (2.36) 


It can be shown that the prediction errors satisfy some orthogonal conditions. These 
conditions are simular to those found in AR modeling problems [Ref. 6: pp. 116-121]. 
We now list the orthogonality conditions for the ARMA formulation in discussion 


without proof as the following: 
EL dk) wA-J) J=0 EL AK) o(kK-) J =0 
El vi{k) x(k-) J=0 Ef Wk) x(kK-) | =0 
ELg.fk —1) HA -/ J =0 (3.37) 
Elig, {A= l) xk —a) J =0 


Elda, (k-—1) y(kK-/)] =0 


Ed, (kh 1) c(h eee 


where (= |(2.2.,7 and) — 1 eee 

In Section B we have obtained a set of normal equations in terms of R,, and the 
ARMA estimation parameters. In this section, we have defined four sets of forward and 
backward estimation parameters and established some orthogonal relationships. In what 
follows, we derive a set of equations which relates the coefficients of the estimiames 
ARMA model with those of the forward estimation parameter matrix C,,. Consider the 
expected value of the forward prediction error, v,,(&) and the data h,,(A). Since the pre- 
diction error is orthogonal to all past samples of data y(k —/), x(k — ‘) but not to 3(k) 


or x(k) as listed in (2.37), the result 1s a matrix which 1s defined as 


T ,y_| si 8 2 8 ; 
EL v,A4) nA®) J= | a i (2.38) 

where 
é=-E[ 40300 | = -E[ thon ]= BF? (2.39) 


is the crosscorrelation between the forward prediction errors of x and y at lag Zero. 

wee pears oval ©) =e, (2.40) 
is the forward prediction error power of x(k) 

é,= E[ 2A) hk) = E| BMY | =e, (2.41) 
is the forward prediction error power of y(k) and 

,=—-E[ 2(Ox(K) )=—-El vate | =e? (2.42) 


is the crosscorrelation between the forward prediction errors of y and x at lag zero. In 


another interpretation, the left hand side of equation (2.38) can be written as, 
E[ v.(k) hk) J= EL Cote (A)” hk) | = Cy Rey (2.43) 


and we have 


cea Was 
Eee, 0 


Ce. Rer= 2.44 
Sf oe vice Ee 0 ited 0 ( ) 
from equation (2.38). 
A similar approach for both backward prediction errors yields the following 
st DES 0 Le 
(2.45) 


te d ed 
D, ; O £,, OE, 


In order to express the coefficients of the ARMA estimation model in terms of the co- 
efficients of the parameter estimation matrix C,, and parameter vectors G,, and D,, , 
consider the combination of the normal equation (2.20) and the parameter estimation 
matrix and vectors. From equation (2.44) the normal equation for ARMA parameter 


estimation may be written as 


Ome Gs Beep, 20 
y y i y xy (2.46) 
1 as, 0 b,, | E;, O £;; 0 
Combining equation (2.46) with (2.20) we obtain 
eA ay ee ie) 0 
UR he eee ee OE. .0 (2.47) 
I as 0 b: Eo D) aia 


ton 


In equation (2.47) multiply row 2 bv = then subtract row 2 from row 3 and equate 


the result to row 1. We then obtain the following relations 











Ese 

as, Fas, 7 as, EX (2.48) 
st 

EN? 
by = = (2.49) 
Ce 

eee 

Do bs, < b;, EX (2.50) 
a 


and, 


(Ex?) 


E,,min = E*, — See (2.51) 
el 


Calculation of the ARMA estimation model coefficients requires knowledge of the est- 
mation parameters a?,, a’,, b%,, bY,, 6, and the values E*,, E¥,, Ex? . Recursive formulas 
calculate the estimation parameters as the AR or MA order of the estimation model in- 
creases by one. Given the parameters of the ARMA estimation model of order (s, 2) 
these formulas calculate the (s+1, 2) or (s, «+ 1) ARMA parameters. For a compre- 
hensive derivation of these recursive formulae see [Ref. 4: pp. 619-621]. The AR-type 
recursive formula for the forward parameter estimation matrix and backward estimation 


paramieter vectors 1s 


T 
Cae =e Cy I, a U; De I, 


ee = L G,, + & D,, ] I, (2.52) 
Dean = Ds; I, os [ U3 oF + Uy G, ; + Us D,, ] I, 
were 
in| 
U = —(2 ee 
ee 
SE 
= =e 
ee ra 
| yr 
U3 ms [ T,1 7% ] Ee (2253) 
aid d d be) 
(ES, time ea) 
= a 
2 d,2 -d 8 
Resse > 6 a 
Us => Us il5 
and the (s+ 1,t) prediction error powers are recursively calulated as follows 
se q] 
Sl e458 l tT, T 
od 
Eo — cee 
Slee 3 2 °4 - 
(2.54) 


d if 


RES 0 gd 
s+1,t 7 ee ag lls Ee 


The matrices J, and I, are of dimension (s+r+2 x s+1+ 3). They are introduced to 
provide symmetry to the matrix algebra and preserve initial condition calculations. We 


design the matrices to perform the following operations 
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L COT le eg rad J I, a [ Dy 0 Wi 0 D542 os Wei ry ] (2) 
[Oy 541 Osta Osarga JI, = 10,» O41) O59 Ossrtt J (2.56) 


The values t, through t, can be obtained through the correlation data and forward and 


backward estimation parameters. We express them mathematically as 


ee ce ners) | 
t73=-ELy(kK-s—1)¢,(k) J Car) 
t4= EL y(k—s—1)d,(k) ] 


The MA-type recursive formula for the forward estimation parameter matrix and back- 
Ward estimation parameter vectors is obtained in a similar manner. The recursive for- 


mula is given by 


ie 
Corel aa Cr I; +N, G. I, 
Doi. = [ Dy, + 1 Gs, i I; (2. 
Gor aaa G,, Hse [ Ng oe + My D, + 5 Cay | I, 


tv 
a 
CO 
an 


where 


h, =— ( ie Les aha. | 


sad 
aye 
ee 
ae 
m=—Ces JES (2.59) 


= (Es 2’, — FE? t's) 
| CES a ar 


Ns ore! i) 15) 


and the (s,t+1) prediction error powers are calculated using the following recursive for- 


mulas 


be = E,, tm, iss eon 
ed =i ' 

ey — Ha + Ith TG (2.60) 
-@) 

P= Get Ot vy }n3 +i elas fits T 


d d od 
Bey = Be, a gh: [So 


where I, and I, are (s + 4+ 2 x s+4+ 3) dimensional matrices which we have designed 


to perform the following operations 


L oo | MG) coo arc j I; = L Oy 2 O54 Oe ee 0 J (2.61) 
L aie ecg) oe kin, coo WO soy J a= L 0 OBER! 0 Os4-2 -2> Ost 549 ] (2.62) 


The values of t’, through t’, are calculated using correlation data in conjunction with 
current forward and backward parameter estimation values. We express these quantities 


mathematically as 


Cr’, ty J=-EL g,fk — 1), f4) | 
3= -EL x(k-0—-1)4,(4)] (2.63) 
t'4= EL x(k—1— 1) ¢,{k) ] 


It 1s interesting to note that the MA-tvpe recursive formula is the complimentary form 
of the AR-type formula and that the two are identical if the variables associated with the 
input signal x(k) and the variables associated with the ouput signal y(4) are interchanged. 
That is, we replace y(k). G,, and g,(k) with —x(A). D,, and d,(A) and vice versa. 
1. Experimental Results 

The ARMA parameter estimation algorithm of [Ref. 4: pp. 619-621] based on 
the recursive formulas of equations (2.52) and (2.58) was implemented using the Fortran 
program found in Appendix A. This program calls subroutines which compute the 
ARMA model parameters as the AR order is increased by one and as the MA order ts 
increased bv one. These subroutines are shown in Appendix B and Appendix C, respec- 
tively. In the main program, an input data sequence of white Gaussian noise 1s passed 
through a known reference model producing an ouput data sequence. We obtain 
autocorrelation and crosscorrelation data from these input and output sequences. The 
correlation data is used to calculate initial values of the error powers for x and y as well 
as t, through t, and 7’, through t’,. Next we obtain estimates of the reference model 
emploving the recursive formulas (2.48) through (2.50), (2.52) and (2.58). Several refer- 
ence models were estimated beginning with a strictly AR process of order s=4 having 


asstSatramsier tuncuon 


| 


QQ) > ———eEee——————————e 
1-022 1406277 SOls2ee 00a 


(2.64) 
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The actual values of the reference model parameters and the ARMA model parameters 


which estimate this reference model are listed below 


ACTUAL SS NGS Eos B) 
AR: 0.2000 0.200583 1 
-0.6200 -0.6207055 
0.1520 Om S27 565 
-0.3016 -).3020376 
MA: 1.0000 1.0001506 


We next consider a second reference model with MA order t=2 and AR order s=3 
having transfer function 


05 = 640s es Sos . 
CO a eee (2.65) 
1—0.202z  —0.252 ~+0.052z 


The true reference model parameters and ARMA model parameter estimates are shown 
to be 


Pelee ESMiviAaleED 
AR: 0.2000 0.1993060 
0.2500 0.2496567 
-0.0500 -).049196]1 
MA: 0.5000 0.5002602 
-0.4000 -0.3997071 
0.8900 0.8894749 


A third example with MA order t= 2 and AR order s = 4 having transfer function 


iets! — ios 
H(z) 


oo (2.66) 
ae G2 ee = 0 15228 0.30162 


was considered for which we obtained the following actual and estimated reference 
model parameter values 


ACTUAL EST ire 
ik 02000 0.2011805 
-0.6200 -0.6225803 
0.1520 0.1534197 
-0.3016 -0.3036823 
MA: 1.0000 0.9997638 
0.2000 0.199§342 
-0.9900 -0.9886852 


We consider as a final example the reference model of AR order and MA order s = 3 and 


t= 3, respectively, with specific transfer function 


f(z) = 
116692 006255 oe 


(2.67) 


The actual and estimated ARMA parameters are 


ACTUAL BSS ATED 
AR: -1.6900 SEGIS 1325 
0.9620 C9075 
-0, 2000 -().2018440 
MA: 0.5000 0.4995653 
-0.9500 “9550507 
1.3300 L533 6707 
-0.9790 -0.9864898 


The above examples demonstrate the validity of the parameter estimation algo- 
rithm of [Ref. 4: pp. 619-621]. Many reference models Were estimated using this algo- 


rithm, including pure MA processes, for which accurate estimates were obtained. 


D. LATTICE STRUCTURE 

In section C we developed expressions for the forward and backward prediction er- 
HoNsmiamely, those oO} equations (2.31), (2.52) and (2.33). From these prediction error 
equations we can design elementary AR, MA and ARMA lattice structures or sections. 
Each elementary section satisfies the orthogonal conditions as listed in equation (2.37). 
From the prediction error recursive formulae, equations (2.31), (2.32) and (2.33), we 
construct the AR-tvpe elementary lattice section as follows. Consider the following data 


set of order (s + 1,r) consisting of input and output data elements, 
ean kyl = [ 3(k)...3(k —s) —x(k)... -x(kK —14+1) —x(k -—d J (2.68) 
he AAI = CAD. K-51) x(k - 1). x(k) 07 (2.69) 


Pere l/ and 17 are the transposes of the matrices I, and I, defined in equations (2.55) 
and (2.56). We obtain a recursive relationship between the forward prediction errors 
v(k) of order (s + 1.1) and order (s,1) by substituting equations (2.52), (2.68) and (2.69) 


in equation (2.31) such that 


18 
V541 (4) - he 1 (4) Cou 
= Thal) 1 EL & Glee ah | (2.70) 
=v, (k) + da, (A — 1) uy, 


The backward prediction error recursions are obtained in a similar manner and the 


AR-tvype error recursions are 


st (A) = i i uy d, {k — 1) 

ae fA) = re (A) + uy me {k — 1) 
10K) = 85 (k) — uy d, (k) 

dk — 1+ [0 8 ] vb? — oy Ben AB) 


Qane 


Gs4y, all ‘} 


where u,=L uy w j and u,=[ w wd. The AR-type elementary lattice inverse section 


based on these error recursions is shown in Figure 5. 
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Figure 5.  AR-tvpe elementary lattice inverse section 


We design the MA-tvpe elementary lattice inverse section in a similar manner using the 


following representation of the data set of order (s.7+ 1) 


eee) i= Ey(A)a(A -— 1). (A — 5) —x(A)... —x(k -— t+ 1) —x(k -—2) J 
hea) (AVI = Lak — 1). (kK -5) 0 x(k - 1). ox(k 1-1] 


and bv substituung equations (2.58) and (2.72) into equation (2.31). we obtain the for- 


ward prediction error recursion as the MA order is increased by one, namely. 


eee, ed (A) + ny ee el) 


a Os : (2.73) 
ern UA) = Vey (k) ky aA f= 1) 


The backward prediction error relationships are obtained in a sinular manner and are 
civen by 
143A) = d, (k) — No g, (K) 


x y T (2.74) 
Sy ar(h) = 8, Ak — ID —[ ny ng vy A) — 1g de aA) 


The MaA-tvpe elementary lattice inverse section based on these error recursions is shown 


In eune Oo 





Fieure 6. MLA-tvpe elementary inverse lattice section 


We now construct the ARMA elementary lattice section from the AR and MA 
prediction error recursions. Assunung that the prediction errors are known for a given 
miearcmorucr (St). the ($7 1.t) prediction errors can be calculated. These prediction errors 
@imercer (s+ 1.t) are then updated as the \iA order increases by one resulting in a pre- 
@renion error of order (s+ 1t+1). We consider the forward predictuon error for X as the 


merwerdcr 1S increased by one. specifically 
es oar x io Was 
ree) tele G de, th — 1) er) 


Seere:?.,,(4) becomes the current value of the fonvard prediction error for x and when 


Lf 


Meemeaiculate the (s+ 1.t+ 1) fonvard prediction error we have from equation (2.73) 
~ a x ° x a yy, 7 
V+ ttl (K) Se Sy (K) a5 My Ssada (K i) ie 6) 


mation (2./G) can be expressed in terms of the (s,t) forward prediction errors of x by 
making appropriate substitutions for »7_,,(k) and g,.,,(A—1). That is, we substitute 


Daeeeand g_,,(A — 1) of equation (2.71) in equation (2.76) to obtain the (s+ 1.t+ 1) 


forward prediction errors for x, namely, 


Vera (A) = VG, (A) — up dd, (kK -— D4 nfl eg, (K-lD-wmd.,(k-1 J (2.77) 


Grouping the terms we obtain 
renee (A) = 0h, (A) — (ul + nf uy) de (A = 1) + af ge (K-12) (2.78) 


The forward prediction error recursion for y of order (s+ 1,t+ 1) 1s obtained in a similar 
manner. We begin with the (s+ 1,t) order update of the prediction error and after it is 


computed, update the MA order. Specifically, we have 


Msi (A) = 154 (K) + Wy de (k= 1) (2.79) 
and from equation (2.73) 
es (k) = ae = m S541, (ha) (2.80) 


Substituting (2.71) for w_,,(A) and g,_,,(k — 1) in equation (2.80) then grouping terms 


we obtain the (s+ 1.t+ 1) fomvand prediction error reciinsiemioi 
ear (A) = 1h (A) + (tq) + re} uy) de (A=W) = at gy (k= 1) (2.81) 


The (s+ 1.t+1) backward prediction errors for x and y are derived in a similar manner 


and arescinen Oy. 


: yy 
Sete = Se (k— 0G Era iayie, Ee eee 


(2.82) 
Asi t+] (A) = Ts (gern U3 ae (is 16 i (k) 


The ARMA elementary lattice inverse section 1s shown in Figure 7 where the coefficients 


are related to the prediction error recursions by the following 


MW) = (ue + wy), wy =n, 3 = (i, +m Us). Ww =a (2.83) 
15 =(n; + Ng 3), iN, = (113 + My Uy); ro ae Wg = iy (2.84) 
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Figure 7. ARMA elementary lattice section 


We see from Figure 7 that each elementary ARMA lattice inverse section contains eight 
eoeiicients. 

From the AR, MA and ARMA elementarv lattice inverse sections, we can obtain 
synthesis lattice structures. These structures provide a means of working with lattice re- 
alizations as linear filters. The resulting AR, MA and ARMA elementary synthesis 


lattice filters are shown in Figures § and 9 respectively. 


Figure 8... Top: AR elementary lattice filter. Bottom: MA elementary lattice section. 
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Figure 9. ARMA elementary lattice filter 


Summarizing. in this chapter we have reviewed the Mullis-Roberts criterion, intro- 
duced the ARMA parameter estimation as a generalized Mullis-Roberts criterion and 
obtained analysis and svnihesis forms of lattice structures. We notice that each ARMA 
elementary lattice section consists of eight reflection parameters and the calculation of 
these parameters requires the autocorrelation and crosscorrelation information as ob- 
tained from the input. output data of the reference model. Also, we obtained a set of 
equations relating the final model estimation parameters and the prediction error model 


parameters which inturn are obtained from the eight lattice parameters. 


to 
tar 


WI. ADAPTIVE LATTICE ALGORITHMI 


A. LEAST MEAN SQUARE ALGORITHM 

The study and design of adaptive filters is known to be a very important part of 
statistical signal processing. Many adaptive algorithms have been developed to support 
the application of adaptive filtering in communications and control [Ref. 12}. An adap- 
uve filter is characterized by the ability of its filter coefficients to adjust (self-optimize) 
automatically and vield an optimum filter design. Two processes occur within an adap- 
uve filter, namely, the adaptation and the filtering processes. During the filtering process 
a desired signal is applied to an adaptive algorithm as a reference for adjusting the filter 
coefficients. Figure 10 shows a block diagram of the adaptive modeling process. Refer- 
ring to Figure 10, let y{A) be the output of the filter at time 4. By comparing the output 
with the desired signal @(A) , an error signal e(A) 1s generated. The adaptive algorithm 
of the filter uses this error signal to generate corrections which are applied to the filter 
coefficients such that an optimum solution is obtained. An optimization technique called 
the method of steepest descent provides an approach to solving this problem. The pro- 


cedure is as follows: 
l. Assign initial values to all filter coefficients. 


2. Using these initial values, compute the gradient vector. whose individual elements 
equal the first derivatives of the mean-squared error with respect to the filter coef- 
ficients. 


3. Compute values for the filter coefficients by changing the initial values in the d1- 
rection opposite that of the gradient. 


4. Return to step 2 and repeat the procedure. 


There is, however, a limitation to this procedure. The steepest descent algorithm requires 
exact measurements of the gradient vector at each iteration which, in practice, 1s not 
possible. Therefore, the gradient vector must be estimated and consequently, errors are 
introduced. An algorithm is required which computes the gradient from the available 
data. The least mean square (LMS) algorithm, developed by Widrow and Hoff, is widely 
used and 1s very convenient to implement in real time hardware [Ref. 13: pp. 96-104]. 
Let }{A) be the output of the filter and d(A) the desired signal at time A as shown in 
Figure 10. We compute the error by taking the difference between these two signals: 


namely, 


d(k) 
DESIRED 
SIGNAL 


AS, y(k) : e(k) 


INPUT OUTPUT . ERROR 
FILTER (+) 
ADAPTIVE 
ALGORITHM 


Figure 10. Adaptive modeling block diagram 


c(k) = d(k) — x(k) (3.1) 


Mivesvalue of the mean-squared error is the expected value of the error squared, 
E{. e(k) | and the gradient vector, V (A) , is the first derivative of the mean-squared er- 


Ree | he cradient Vector 18 given by 


CEN ‘i {ry — io a. 5 
V(A) = TTS E| e*(k) | =2 (A) En e(k) (3.2) 





Where w(k) is the time dependent filter coefficient vector. The recursion for the filter 
coeflicient which changes the old value in the direction opposite to that of the gradient 


ieecnen given by, 


wk)=wk—-N++ al -V(A)] 
yas | ne . 
= w(k — 1) — yw e(k) Ue e(k) 
Where w(A) is the filter coefficient vector estimate at the A“ iteration, w(k — 1) is the past 


filter coefficient vector estimate, y is the convergence (gain) constant, e(k) is the error 


~ 


° je 4 C 
signal at the k” iterauon and 
cw(k) 


tation of this algorithm proceeds as follows: 


e(k) is the instantaneous gradient. The implemen- 


1. Assign initial values to the filter coefficients. 
2. Compute the vwaluewer the error sicital efk 


3. Calculate the updated estimate of the filter coefficients using the instantaneous 
gradient. 


4. Increment the time index bY one and Tetum to sitcom 


Convergence properties of the LMS algorithm are well documented within the literature. 
The choice of a gain constant yp is arbitrary however, theoretical bounds have been de- 
rived for np, given by [Ref. 14: pp. 101-106], 


| 


a 
O<u <2 3.4 
a Trae (3.4) 


“max 
where /,,,, 18 the maximum eigenvalue of the input autocorrelation matrix, R,,. and 


where Tr [ R,, | is the trace of the matrix R,,. 


B. DERIVATION OF THE ADAPTIVE LATTICE ALGORITHM 

The adaptive lattice algorithm developed in this thesis uses concepts of the LMS 
algorithm discussed in section A and applies them to the ARMA digital lattice filter 
proposed in Chapter II. Consider the ARMA digital lattice filter of Figure 11, which 
consists of two cascaded elementary lattice sections. The filter coefficients (weights) are 
defined such that 1” represents the i” lattice coefficient at stage + of the lattice structure. 
In this figure we have a two stage lattice and there are eight coeflicients per elementary 


lattice section. The output, s(A), of the lattice filter can be determined from 


or) 
Cay 
ee” 


B(k) = 6 (k) + we e5,(k — 1) — wy eg (k — 3) (3. 


Forward errors at a given stage m of the lattice filter are defined as, 


ez (k) = er (k) + wy Cy (kK-l)—w) ee (k-1) 


| | _ 3.6 
gf, (k) = Gf. (H) + ut ef (k= 1) = ws (k= 1) Ce 


and the backward errors for any given stage mm are, 


cy (ile ey (A—1)+ we" er (k) — we cp (k) 
e (k) = crm (k—-l)—we oa (kK) + wy" cp (A) 





—_ 


Two stage ARNIA lattice digital filter. 


Ficure 11. 


Where, in Figure 1]. m= 1,2 and ex =x(k) and e} =y(k). The terminal condition is 
ey (kK) = bef (kK). m=2 . To begin with, let 5, equal unity. The initial conditions are 
ey. (kK — 1) =0 and e; (A — 1) = 0. As with the LMS algorithm. we form an error between 


a desired signal d(k) and the output signal +(A) such that, 
e(k) = d(k) — ¥(4) (3.8) 


The instantaneous gradient according to eq (3.2) is then, 
Nae Cm ae 3 
Vik) = 2 ek) Bgy Late) - A) J (3.9) 


Since the desired signal, d(k), is not a function of the filter coefficients, equation (3.9) 


reduces to, 





V(k) = 2 e(k) ay [ -¥(4) J (3.10) 


where the quantity [ —y(k) ] is referred to as the gradient estimator. This gra- 


C 
eWw(k) 
dient estimator must be computed for each filter coefficient within the lattice structure. 


The filter coefficients are then updated using the respective gradient estimators. That 1s, 


we need to compute. 


~ 


Vik‘) = a ? (kh) fori=1.2,...,8 and j= 1,2, 20a (3.11) 





CAs 


where M is the number of stages in the ARMA lattice filter. From equation (3.5), the 


gradient estimates are given by 
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Let w(K) represent the partial derivatives of the output }(k) with respect to the filter co- 
efficients and $(k) represent the partial derivatives of the errors with respect to the filter 


coefficients. Using this representation we can re-write equation (3.12) as follows, 
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bij (kK) = bp i(k) + O54 ey, (A — 1) + 9g 6517 (K — 1) — 634 2, (K— D) 


=; Dy, ij (A i 1) ss 


— 
Lad 
Nee” 


where olf, d!/, are kronecker deltas whose value is one if and only if 
i=4 and j=1 or i=3 and j= 1 respectively. We compute y,,(A) , by obtaining re- 
cursive relations which calculate the partial derivative of the forward and backward er- 
rors with respect to the filter coefficients. These relations are obtained from equations 
(3.6) and (3.7) by taking the partial derivatives with respect to the filter coefficients, 


namely, 


OF i(k) = OF aj Ck )+ 55) ex (kK-D+ wy by i (K-D- 5TH ep (ee, 
—Wy Ps, me I) 


BF i(k) = OF a (k) + ONT eg (kK V4 wl! OF (K-D— OS @f (k= 1) 


m+] 
me Oy, ij (a) 


6, (M=o5 (kK-N + 85/ ef (A) + ws OF is (A) — 867 of, (4) 
We oF, ij (A) 


Pp, y (k) - Py ij (k 7” 1) ia 577 er 1 (k 6 <i 1 Dy yi (K) +r ee Cp ; (K) 
+g" op, a (x) 


These recursive relations possess a lattice structure similar to that of Figure 11. with 
delta components injected at the summation nodes. Thev mav also be simplified bv ex- 
anuning individual terms. Consider the general equation for the forward error in x, re- 


peated here for continuity. 


er (k)= er (kK) + wy eee (kK-—1) — wy Fw (k—1) (8 1p) 


The partial derivative with respect to each filter coefficient is expressed as, 





malpee,4 A 
Cer (k) Cer (Kk) Brees - Z Ge, (A — 1) 
a j ie Chik I) + 5 —— 
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Since the partial derivative is taken with respect to the current filter coefficient w! at time 
A, the partial derivatives involving delay terms 1.e., (k — 1). are set to zero. This result 
follows from the realistic assumption that e; (A — 1) is a function of w!(k — 1) but not 
of wi (A) . Also note that w’ (A) is a function of w(A — 1) but not viee versa. With these 


simplifications we reduce the equations of (3.14) to, 


OF (A) = O71 (A) +6 Oo) ey. b= Wes 11 (kK 1) 
OY (KV = OP_ ij (A) + OY! of (k- 1) ged gp (k- 1) 


sid ox (Sl 
bs, 1p (K) = 05 557 Ge Wy) oe Ws oF 1 (k) — 667 {ef (A) — we" oy, ij (A) 
bp, i} (A) = — 67) & ef) (A) — Wy" oF, tj (ee Ont er (A) + We op ij (x) 
and the gradient estimator 1s. 
Wiz (A) = OF 1 (A) + 524 @6, (A — 1) — 534 €2, (K-10) (3.18) 


Although these are valid recursive relations, they are difficult to implement in a lattice 
algorithm. The ultimate goal is the requirement to easily compute w,,(k) from the 
available data. From eq (3.18) it is evident that y,, (A) depends on @} ,, (A) which in turn 
(A) but not of @; ,,(A) or @} ,, (kK). Therefore the 


three equations necessary to compute the gradient estimator are, 


requires knowledge of py ,, (A) and @; 


fmt Omtf 


Wig (A) = bp aj (A) + 544 ep, (A — 1) — 637 5, (A — 1) 

OF i (K)=O7 (A) + 677 eo (A-N- 67 eg (A- 1) (5.19) 

Or (A= Op i(k) t+ OE of (k-1)- sed ep (k- 1) 
These equations are dependent on the filter coefficients w/, 1=1.2.3.4 and 
j= 1,2,....3f, thereby reducing by one-half the number of computations required for 
ob; .,(k) and @ ,,(A). A recursive relation is desired for p,,(k) which does not involve 
delta functions. Consider the four stage lattice filter with terminal condition 4, equal to 
unity such that },,;(k) = $3,,,,;(k) . The procedure for computing w,,(k) using the 
equations of (3.19) 1s as follows: 


1. Calculate $; ,,(k) with mz equal to one and letting 7 and j range from one to four. 
Repeat for m= 2,3,8. 


2. Using the terminal condition and expression for @} ,,(k), calculate @+ ,, (4). 


This requires solving 88 equations, however, the result is a verv simple recursive formula 


for w,,(k) , namely, 


Wij(=-e (k-1) i= 1.3 — 
: aie 

Wij(K=ek (kK-1) i= 2,4 ) 
The lattice coefficients are calculated bv substituting the recursive formula of eq (3.20) 


for the gradient estimator in equation (3.3). The adaptive coefficient update equation Is, 
w7 (k) = 07 (K — 1) — 2 eK) Wr (4) (3.21) 


Since the gradient estimator, and therefore the gradient, is the same for w , i= 1,3 and 
wi, i= 2,4, it follows that wi = ws and w! = wi. The number of filter coefficients re- 
mimired' to Update the lattice filter is reduced to two, 1e., wi and wi. Furthermore, from 
the symmetry of the lattice structure, the following equalities between filter coefficients 


are assumed, 


wm = WwW 

i — ae ‘ 
j (Sea) 
| a E 

Ws 4 
J “| 

Wz = Wy 


Incorporating these equalities with those derived by the gradient estimator produces the 


elementary ARMA lattice section of Figure 12 where. 


ne = Wi, = We = y= i 
4 55 
[oa 
“7 = Wr = We = he — ik , 


To prove that these coefficient reductions are valid, a computer generated solution using 
the Fortran program of Appendix D was compared to hand analysis of a second order 
transfer function and lattice filter. The output of the ARMA digital lattice filter was first 
put into difference equation form and then compared to the known transfer function. 
From this comparison, lattice coefficients were computed. Details of this analvsis are as 
follows. Consider a two stage ARMA digital lattice filter comprised of the reduced ele- 


mentarv section shown in Figure 13 and a transfer function of the form. 


bo + b, 7a = b, as 


—| —2 
1+ a, z aie toe 
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Figure 12. Simplified elementary ARMA lattice section. 


The output of the lattice filter can be written in difference equation form by carrving out 
the following steps: (i) start with the output of the lattice filter equation (3.5), (11) sub- 
stitute expressions for the forward and backward errors, equations (3.6) and (3.7). re- 
spectively. into equation (3.5) and (il) carrv out the algebra. A detailed derivation of this 
difference equation is given in Appendix E. The difference equation in its final form 1s 


given by, 


(hk) = 2x(4) + 2 5 - "y Wy ~ 3 “5 )x(A — 1) +2 re x(k — 2) 
a 


—2( Wy -- ve ws “+ wy We yk — 1) — 2 wy ¥(k — 2) 


which can be written in the transier 1unction om 4s. 
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f(z) 


] as2 eee ==] 2_- 
1+2( 4, + wy) wy) + wy wy )z 5 +2452 


~ —_ a 
P+ 20, + wo wy tow, wp )z +2) 2 


2 


= (3.26) 


~~» 
No 


Comparing the lattice filter transfer function, //{z) with the known filter transfer func- 


tion Hz), produces the following relationships between filter coefficients, 


» 


ee ne: 
b, =2( Ww, + wy, wy + wy wy) 

oes epee 
a= 2( wy) + Wy, wy + Wy, Wy) 297 
—_— (G27) 
Ae AD. 
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Solving for #3 and 7? in terms of the known transfer function coeflicients 6, and a, and 


then substituting these results into the expressions for b, and a,, respectively, vields the 





following. 
by = aw, + (2+) 0 256 
ti pes J sl (5.28) 
a, =(2+a,), + by Ww 
rein matrix form, 
a yy Ta Ww 
. “Gis a b : . C4 
- 7 ay 2 2 
and solving for the lattice coefficients, we have 
Wy 1 b, — (2+ b) b, (3.30) 
— a a =) 
5 ay b, _ (2 oF b,) (2 Ir ay) = (2 ai A>) Qy ay 
and 
2 
Dale: 2 
a Ga (r3)) 
WwW) = 5 


Now that a method of converting between lattice and transfer function coefficients for 


a second order system has been established, we consider the specific transfer function 


l= 082° 2 ieee 


Hz) = 
f 1— 0892" = 0 257m 


~~ 
{od 
or) 
tJ 
—Y 


where 


by =1.0 b,=-0.80 b,=1.78 
-~0.89 a, =0.25 


4 
From equations (3.30) and (3.31) the lattice coefficients are calculated as, 
we = 0.125, wy = 0.890, w, = —0.240719, wy = —0.195719 


Values for the steady-state lattice coefficients were computed using the Fortran program 
in Appendix D and are shown below. Convergence aspects of both the lattice coeffi- 


cients and output error are shown 1n Figure 13. 
wy = 0.124982, wi = 0.890003, wi =—0.240710, wy, =—0.19571] 


these results confirm the validity of the derived adaptive lattice algorithm and the design 
of a new elementarv lattice section shown in Figure 12. 

The current adaptive lattice algorithm assumes that the terminal condition is unity. 
This is generally not the case in practice. We now extend this adaptive algorithm to the 
more general case where the terminal condition is an arbitrary constant. The recursive 
relation Which updates 4 is similar to those which update the other lattice filter coeffi- 


cients. The update equation for }, is given by 











bok) = blk — | joes 338 
of )= o( ar )— pu ef ) 6ho(k) (3.35) 
| . Ce(k) . : 
The gradient estimator bak) is calculated using equations (3.5), (3.8) and the fact that 
the desired signal d/k/ 1s not dependent on 4. The gradient estimator for b) is written 
as, 
A ae ae a | Ane ~ 1 
Cy(k Ce, (k) Cw Cen (coum) Cw 
scl) Meee _ é, (ae eee ———  - —— e,(k — 1) 
0 bo C bo Cbg 0 Cho Cho 0 


(3.34) 


1 
—— sh 3 


de, (k —]) 
Cho 


Since the partial! derivative is taken with respect to 5, at time A, tis senUeessme: 
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Figure 13. Top: Lattice coefficients. Bottom: Output error. 
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Eb s aa, (3.339 
Simularly, 
ep (Kyser (A) + wy eg (K-11) — wy"! ef (k-1) (3.36) 
and the partial derivative with respect to 6, 1s 
la Cae Oo (3.37) 
Cho Cho 
The terminal condition is, 
er (k) = by ep (k) (3.38) 


Taking the partial derivative of equation (3.38) with respect to 4 yields e?, (A), and the 


recursive equation to update b, (k) becomes, 
bo(k) = bok — 1) — we(k) ef, (A) (3.39) 


The gradient estimators for the lattice filter coefficients are scaled by the arbitrary con- 
(kK) = b, @} ,,(A) and 6 1s propagated 


through the calculations. The gradient estimators become, 


stant 6), since at the terminal condition @},, 


Wij(K)=—boeh_(kK-1) i= 13 | 
(3.40) 

W;; (A) = by cs = Ak = ieee 

To test this more general adaptive lattice algorithm the output of a known transfer 

function with & equal to 0.5 was compared to the output of the ARMA digital lattice 


filter. The second order transfer function used was, 


05 =e ee nc ee ; 


Hz) = (3.41) 
: 0/90 es 


The computer generated steady state lattice coefficients are given below and convergence 


aspects shown in Figure 14. 


by = 0.499946, wy = —0.120428. wy = 0.593237, w, =—0.447423, wy = 0.166320 
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Figure 14. Top: Lattice coefficients. Bottom: Output error. 


ay 


In order to maintain the same adaptive time constant and misadjustment at each 
stage in the lattice, the convergence constant is normalized by the power level at each 


stage [Ref. 15]. Therefore, we can write equations (3.21) and (3.39) as, 


ei (3.42) 
by(k) = Belk — IJ — e(k) 6 (4) 
y(A) 


wl (k) ow! (k= De Ce 
; (A) = wy; ( ) 2h) e(k) W;;(k) 


Where yw 1s the convergence constant and a?(k) and y?(K) are estimates of the power at 


the j" stage for w’ and 4, respectively and computed as follows: 


oj (k) = p oj (k— 1) + (1 p) Wiz (K) 


3.43 
y (A)=py (kK-1)+(1-p)] && (4) | i 


Writing equation (3.42) using the notation adopted for the reduced elementary ARMA 


lattice section we obtain 


r(k) =r (k- 1) a e(k) ef (k= 1) 


a 
Ll ae 
kj (k) =k (K- D-—-—— ek) (k-1 4 
(3 (A) i A ) G7 (k) e( ) pe ) (3.44) 


by (k) = By (k — 1) —— ek) ef, (A) 
y” (A) 
In the above equations p is a Weighting parameter, 0 < p < 1, which distributes the 
amount of weight given the past power level or current sample. Normalized convergence 
constants are used in all examples of this thesis. The adaptive lattice algorithm 1s sum- 
marized in Table |. 

In summary, we have derived an adaptive algorithm based on the LMS theory of 
adaptive coefficient computation. This new adaptive algorithm easily updates the lattice 
coefficients by using available data. The original requirement to update eight coefficients 
of an elementary ARMA lattice section was reduced to updating onlv two coefficients 
and still being able to describe the lattice. The algorithm is general in that it apphes to 
svstems Whose terminal condition is an arbitrary constant. The validity of this algonthm 
was demonstrated through comparisons between hand analysis and computer simu- 


lation. In the next chapter, we further demonstrate the convergence of this algorithm. 
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Table 1. SUMMARY OF ADAPTIVE LATTICE ALGORITHM 


With given initial conditions, e7 = x(k), ef =y{k), ef (A — 1) = ef, (A — 1) =0 and all 
lattice coefficients zero, 


Step I: for k=1,m compute 


(k) + w5" ey (A—1)—w, ey (A — 1) 
(kK) + 1 ey (K-—1)- wit! eh (k — 1) 


m—1 


with output e} (k} 


Sicep 2: for k= 1,m compute 


Ci, Co — ca (A — 1) + we ep (k) — we er (k) 
ey (‘ie coe (K-11) — wy op AGS an Wg: or (A) 


elep >: Update coeilicients 


1 (K) = 1) (k — 1) —-— eh) ef (K-1) 
oj (k Z 
Lu 
k.(k) =kj(k--— 
ai 
U 
bo(k) = bo(k — 1) — — eA) e7-_ (A) 
yk) 


e(k) apa, (A — 1) 





meg Repeat for next iteration i.e. return to step |. 
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IV. EXPERIMENTAL RESCUES 


The adaptive lattice algorithm derived in Chapter III 1s now computer simulated to 
study its convergence performance. The system identification mode of adaptive filtering 
is considered for this purpose. Figure 15 shows a general system identification config- 
uration. The systems considered are time-invariant and linear. Notice that We apply the 
same Input, white noise in general, to both the reference system and the adaptive lattice 
filter which 1s modeling the system. The criterion in this configuration 1s to minimize the 
miean-squared error between the system and filter outputs. Thus, in this context, the 
adaptive algorithm continuously updates the lattice filter parameters in order to mini- 
nuze the mean-squared error. 

The adaptive algorithm is realized as summarized in Table 1. As We mentioned in 
Chapter III, the two important parameters of the algorithm are the adaptation constant 
u and the weighting constant p. In what follows, we shall! consider convergence studies 
of both second and third order reference svstems (fixed filter transfer functions). Con- 


sider the following reference system with transfer function, 


[- 20225 = Oesee 


Hz) = 
g = iee 4.10 Cor 


This system has complex poles and simple zeros located at z = (0./+/ 0.6) and 
z = 0.5, —0.7, respectively. Using a convergence constant » = 0.01 and power level 
Weighting factor p = 0.45, the adaptive ARMA digital lattice filter which models the 


above systein has the following steady-state lattice parameters. 


terminal condition by = 0.999202 


lattice coefficients r; = 0.352580 
r7 = —0.174355 
ki = 0.448228 


k? = 0.425060 


il 


Convergence properties of the lattice coefficients and error are shown In Figure 16. The 
mean-squared error was minimized after approximately 1700 iterations at which time the 


lattice coefficients reached their steady-state values. When the value of the convergence 


reference 


system 





Figure 15. Block diagram of system indentification modeling 


constant un was modestly increased, convergence degraded rapidly. Also, when the 
Weighting factor p was increased, convergence deteriorated quickly. From this, we con- 
clude that the convergence constant 1s the more sensitive input parameter. 


Let us consider another second order dynamic system with transfer function, 


_ OS —02e POEs = = 


fi ke 
i en oon 


iemesvstem has complex poles at z = (0.5+/0.8) and complex zeros located at 
z = (0.2 +/0.7). Using a convergence constant 4 =0.005 and power level weighting 
factor p =0.97, the adaptive ARMA digital lattice filter which models this svstem has 


steady-state lattice parameters as follows, 


terminal condition by) = 0.500027 


0.104654 
r? = 0.296658 
ki = —0.429232 
ki = 0.627718 


— ! 
lattice coefficients 7, 
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Figure 16. Second order ARMIA lattice filter, terminal condition unity. 
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Convergence properties of this adaptive filter are shown in Figure 17. In this example, 
convergence was obtained after approximately 2500 iterations. Again, by changing the 
values of uw and p slightly, covergence deteriorated with the convergence constant yu being 
the more sensitive parameter. 


Next we consider a third order reference system with known transfer function, 


Az) = - a z 
1— 1.6927 +0.9622-° — 0.22 


The adaptive ARMA digital lattice filter which describes this system has the following 


steady-state lattice parameters, 


terminal condition 6) = 0.499970 


lattice coefficients rj = —0.328447 
r? = 0.399472 
r= —0.652706 

—e-0.821738 

= —0.091111 


= —0.133335 


These parameters were obtained using a convergence constant un = 0.015 and power level 
weighting factor p=0.9. Convergence properties of this adaptive filter are shown in 
Figure 18. Steady-state values for the lattice coefficients were obtained after approxi- 
matelv 7100 iterations. It is reasonable to assume that a third order svstem will converge 
more slowly than a second order svstem. The number of iterations required for this third 
order system to converge is consistant with convergence rates of other adaptive algo- 
rithms which model third order systems [{Ref. 16]. The input parameter uw was again 
found to be the more sensitive parameter. 

In all the previous examples, the values of uw and p may or may not be optimum 
values. That is, an exhaustive search of all combinations of u and p was not performed 
to demonstrate convergence of the algorithm. Nevertheless, a number of different wavs 
of realizing the value of the convergence constant yw have been reported in the literature. 
In one method, Mikhael et. al. [Ref. 17] have obtained a variable uw by using a self opti- 
nuzing technique. In this method, yw is calculated from the input data as an iteration 
process and is individually determined for each filter parameter. In another method yp 1s 


chosen by using a variable step LMS technique [Ref. 18]. where the range of yw 1s 
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Fiecure 17. Second order ARMA lattice filter, terminal condition of 0.5. 
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Figure 18. Third order ARMIA lattice filter, terminal condition of 0.5. 
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specified by u,,,, and yy,, Which are within the bounds described by equation (3.4) of 
Chapter III. These techniques of choosing yu during the adaptation process have been 
shown to inprove filter convergence. They, however, require additional computations 
to achieve this faster convergence. When a combination of the parameters yp and p was 
obtained which vielded convergence, these values were chosen for examples. Besides the 
examples reported, simulation studies have been carried out for several other cases. In 
all cases, however, definite convergence of the algorithm has been observed. 

In summary, we have demonstrated through computer simulation that the derived 
adaptuve lattice algorithm 1s suited for svstem identification modeling. Furthermore, we 
have shown that there 1s flexibility in choosing the values of the convergence constant 
uw and weighting factor p . Some techniques for selecting (computing) the value of the 
convergence constant have been introduced. These methods improve convergence at the 


cost of additional computations. 


A. CONCLUSIONS 

In this thesis we have demonstrated that the ARMA parameter estimation algo- 
rithm proposed in [Ref. 4: pp. 619-621] is a valid method for obtaining approximations 
to reference models. Furthermore, the criterion used to derive the algorithm 1s a gener- 
alized form of the Mullis-Roberts criterion for least squares modeling. The AR and MA 
parameters of the ARMA model can be updated independently as their respective orders 
increase by one. From the recursive prediction error formulas, an ARMA digital lattice 
filter was designed with arbitray AR and MA orders. 

For the ARMA digital lattice filter, we derived an adaptive lattice algorithm. This 
algorithm was based on the least mean square method of optinuzing coefficients. The 
derived adaptive lattice algorithm can easily compute the values of the lattice coefficients 
from available data. The algorithm simphfied the number of coeflicents required to be 
updated from erght coeflicrents per elementary lattice section to only two such that the 
filter can be completely described. This savings in computational effort makes the al- 
gorithm attractive for identification of unknown systems since many systems require an 
ARMA model for parsimonious modeling. 

Convergence of the adaptive lattice algorithm was demonstrated with several exam- 
ples in Chapter III]. The number of iterations required before convergence varied greatly 
between second and third order models as well as within second order models. Optrmum 
convergence rates were not sought after as much as proving the convergence of the al- 


gorithm. Rapid convergence rates were demonstrated in Chapter II for a second order 
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system upon completion of an extensive search for the optimun values of the conver- 
gence constant and power level weighting factor. 

Although a method of converting between direct form and lattice realizations for 
second order ARMA filters was developed, a general algorithm to perform this trans- 
formation given the ARMA lattice filter design was not obtained. When solving for a 
transformation between filter realizations of third order the solution is hindered by 
nonhnearities. 

The objectives of the thesis were sucsessfully accomplished. Some suggestions for 
future work include the following: (1) extensive theoretical analvsis for determining op- 
timum values for p , (11) derivation of a generalized algorithm which converts any given 
ARMA transfer function into a set of lattice parameters, (in) development of theoretical 
convergence models for the ARMA adaptive lattice algorithm and analvsis of these 
models and (1v) apphcation of the adaptive lattice filter, both analysis and synthesis 
forms, in modeling such practical signals as speech. ARMA lattice filter modeling has 
considerable application potential because of its very accurate modeling of nearly any 


signal or system of interest. 
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APPENDIX A. MAIN PROGRAM TO ESTIMATE ARMA PARAMETERS 


Can) Glue, Gia Ge) Cyc) OC) C2) Gl C).C) CoC) Cy Cy) C2 CCG) 0) 0 ea GC) C0 6) Cae) C0 C3 - 02 0) a eG) Ge) © a Cae 


THIS PROGRAM COMPUTES THE ARMA PARAMETERS AS THE AR OR MA ORDER 
OF AN ARMA MODEL INCREASES BY ONE. IT USES THE ARMA PARAMETER 
ESTIMATION ALGORITHM PROPOSED BY MIYANAGA, NAGAI AND MIKI. 


VARIABLE DEFINITIONS 


VN - INPUT VECTOR CONTAINING DATA GENERATED AS WHITE GAUSSIAN 
NOISE 

4 - OQOUPUT VECTOR OF DIFFERENCE EQUATION WITH VN Ase aNEGs 

RX - AUTOCORRELATION DATA OF INPUT VN 

Jsio - AUTOCORRELATION DATA OF OUTPUT Y 

RXY - CROSSCORRELATION DATA OF INPUT AND OUTPUT 

RYX - CROSSCORRELATION DATA OF OUTPUT AND INPUT 


NDATA - NUMBER OF INPUT DATA POINTS 
KDATA - NUMBER OF INITIAL DATA POINTS TO DISREGARD 


X1 - INPUT DATA VECTOR AFTER DISCARDING KDATA POINTS 
pas - OUTPUT DATA VECTOR AFTER DISCARDING KDATA POINTS 
A - VECTOR CONTAINING CALCULATED AR PARAMETERS 
B - VECTOR COINTAINING CALCULATED MA PARAMETERS 
TXA - VECTOR CONTAINING COEFFICIENTS FOR FORWARD PREDICTIONS@ 
ENP UT 
TYA - VECTOR CONTAINING COEFFICIENTS FOR FORWARD PREDICTION OF 
OUPUT 
TXB - VECTOR CONTAINING COEFFICIENTS FOR FORWARD PREDICTIONS 
IS elONe 
TS - VECTOR CONTAINING COEFFICIENTS FOR FORWARD PREDICTION OF 
OUDPUE 
GA - VECTOR CONTAINING COEFFICIENTS FOR BACKWARD PREDICTION OF 
INPUT 
GB - VECTOR CONTAINING COEFFICIENTS FOR BACKWARD PREDICTION OF 
TNR 
ZTA - VECTOR CONTAINING AR COEFFICIENTS FOR BACKWARD PREDICTION 
OF OUTPUT VALUES 
ZTB - VECTOR CONTAINING MA COEFFICIENTS FOR BACKWARD PREDICTION 
OF OUTPUT VALUES. 
NI = LENGE OF DATA VECTORS! 
NK - LENGTH OF DATA VECTOR X1 MINUS ONE USED TO START 
CORRELATION COMPUTATIONS. 
NS - DESIRED AR ORDER SOr ART ASHOM iit 
NT - DESIRED MA ORDER OF ARMA MODEL 
KS - CURRENT AR ORDER OF UPDATE 
ay - CURRENT MA ORDER OF UPDATE 
VX - EXPECTED VALUE OF PREDICTION ERROR FOR INPUT SQUARED. 
VY - EXPECTED VALUE OF PREDICTION ERROR FOR OUTPUT SQUARED. 


VXY - EXPECTED VALUE OF PRODUCT OF PREDICTION ERROR FORKS Eaa 


ts 
= 


ae Co Cl ten Co C) 


>) 


ap Ties Sp OS See a SP a 


10 


LS 


20 


ANDSOUTPUT. 


VG = Ger EOlED VALUE OF BACKWARD PREDICTION ERROR OF INPUT 
SQUARED. 

VZ - EXPECTED VALUE OF BACKWARD PREDICTION ERROR OF OUTPUT 
SQUARED. 


VGZ - EXPECTED VALUE OF PRODUCT BETWEEN BACKWARD PREDICTION 
ERRORS OF INPUT AND OUTPUT. 


DEE NomONS VNC= 57006) 21(-5: 2006) ,RACO: 2000) ,RY( 07 2000) ,RXY( 0: 2000) 
Divi NomONmhieCOr2 000 je x1( 2000), Y1IC 2000) ,XC 12) ,A( 0221), BC 0: 21) 

DIME Swen EACO: Zi ye Ty Ar 0:21), UxXB( 0: 21) ,TYBCO: 21) ,GACO0: 21) 
DEMENSION GEG 2 ZTACO: 21) ,ZTBCO: 21) 


INPUT DATA INFORMATION 


Wie (One) 

FORMAT (/' ENTER THE NUMBER OF DATA POINTS: ') 

READ (6,%*) NDATA 

Wht ti eCon2) 

FORMAT (/' ENTER NUMBER OF INITIAL DATA POINTS TO DISREGARD: ' ) 
READ (6,*) KDATA 


+ eps 


tesevevodevodesevevedesoucsevevovedesedcsevesevedesosescsevesescsedevedesleseseskseslesicvlesicsvvacvevercsevescvesesvesevescdcseveveveveve 


DO 10 L=-5,NDATA 
VN(L)=0 

Y(L) =0 
CONTINUE 

DO 15 L=0,20 
A(L)=0 
B(L)=0 
TXA(L)=0 
TXB(L)=0 
TYA(L)=0 
TYB(L)=0 
GA(L) =0 
GB(L) =0 
ZTA(L)=0 
ZTB(L)=0 
CONTINUE 

DO 20 L=0,1999 
RX(L)=0 
RY(L)=0 
RXY(L)=0 
RYX(L)=0 
CONTINUE 

DO 25 L=1,12 


C2 Go Coach 


LS) 
7a) 


oO 


Ee Ms SP FT Sap a Ga Sk Ss 2 Ba 


Ca Coe 6) ho 6) Cy) 


eee coe) Ga CuGeO) Ct) 0) C) 06) Co 


X(L)=0 
CONTINUE 


Sevedovedodedsuodevededosedodedevesededevlededvvcdedesevevededevedcdedededededededede deve devedcdedevevcdvdevededesccstccssesestesesesteoe 
GENERATE WHITE GAUSSIAN NOISE INPUT 


ISIZE=NDATA 
TX=152255 
ISORT=0 
MUL=2 
DO 30 K=1,ISIZE 
CALL SRND(IX,X,12,MUL, ISORT) 
XT=-6. 0 
DO's Sie ie 
XT=XT+X(T) 
VNCK)=XT 
CONTINUE 


seacvencicassvesedevevevcseveslcveveskvsevlesdesicvevestvesevescsescviesesievesescsevevcvevcdedcvlevesess deslcvlecslevcvcuesvovdesevleveseveveveve 


COMPUTE OUTPUT OF REFERENCE MODEL FILTER AND DISREGARD SPECIFIED 
NUMBER OF DATA POINTS 


DO 40 L=1,NDATA 

Y(L)=VN(L)41. 6*°Y(L-1)-0. 95*Y(L-2) 

Y(L)=VN(L)+0. 2*Y(L-1)-0. 62*Y(L-2)+0. 152*Y(L-3)-0. 3016*Y(L-4) 

Y(L)=VN(L)-0. 2*VN(L-1)+0. 62**VN(L-2)-0. 152*VN( L-3)+0. 3016**VN(L-4) 

Y(L)=VN(L)+0. 2**VN(L-1)-0. 99*VN(L-2)+0. 2*°¥(L-1)-0. 62*Y(L-2)+0. 152*Y 
&(L-3)-0. 3016*¥(L-4) 

Y(L)=VN(L) -1. 6*VN(L-1)41. 45*VN(L-2)+1. 2*Y(L-1)-0. 72*Y(L-2) 

Y(L)=VN(L)+0. 2*VN(L-1)-0. 35°*VN( L-2)41. 4*Y(L-1)-0. 85*Y(L-2) 

Y(L)=0. 5°*VN(L)-0. 2*VN(L-1)+0. 445°°VN(L-2)+Y(L-1)-0. 94*¥(L-2) 


Y(L)=VN(L)-2. 7*VN(L-1)+3. 21*VN(L-2)-1. 595*VN(L-3)+1. 95*Y(L-1)-1. 62 
&*Y(L-2)+0. 54*°Y(L-3) 

Y(L)=VN(L)-1. 0**VN(L-1)+0. 89*VN(L-2)+0. 40*Y(L-1)-0. 2121*Y(L-2)-0. 20 
&894*Y(L-3)-1. 810373*Y(L-4) 

Y(L)=0. 5*VN(L)-0. 95*VN(L-1)+1. 33*VN(L-2)-0. 979°*VN(L-3)41. 69°°Y(L-1) 
&-0. 962*Y(L-2)+0. 20*Y(L-3) 

Y(L)=0. 5°*VN(L)-0. 4**VN(L-1)+0. 89*VN(L-2)+1. 69°*Y(L-1)-0. 962*Y(L-2)+0 
&. 2*Y¥(L-3) 

Y(L)=0. 5°*VN(L)-0. 4*VN(L-1)+0. 89° VN(L-2)+0. 2*¥(L-1)+0. 25*Y(L-2)-0.0 
&5°*Y(L-3) 

Y(L)=0. 5*VN(L)-0. 4*VN(L-1)+0. 89*VN(L-2)+0. 89*Y(L-1)-0. 25*Y(L-2) 
Y(L)=. 0154%*VN(L)+. 0642*VN( L-1)+0. 0642*VN(L-2)+0. 0154*VN(L-3)+1. 99% 
&Y(L-1)-1. 57*Y(L-2)+0. 4583*Y(L-3) 

Y(L)=0. 5*VN(L)+0. 256*VN( L-1)+0. 1234*VN( L-2)+0. 0987*VN(L-3) 
CONTINUE 


{fr 
ty 


200 


Garo) OC) 3° C27 oO) 


LJ=KDATA+1 
DO 45 L=LJ,NDATA 
LK=L-KDATA 
X1CLK)=VN(L) 
Y1CLK)=Y(L) 
CONTINUE 
WRE PEC 37 7) Gr mis, K=200, 1800) 
FORMAT SPA Oro) 


seovcsesicsvedescskvdevedeslescvescsesdedkedeskukscsevevedk seskleskvse dese se seskotakeakskeseveveslesestckvvkvdesese obsess seve sk seskvevese eve 
COMPUTE AUTO-CORRELATION AND CROSS-CORRELATION TERMS 


NI=NDATA-KDATA 

NK=NI-1 

GAMMECORREL (NI,50,X1,Y1,RX,RY,RXY,RYX,NK) 
Dees oml—o- 10 

Neti 200) RXCL) RYCL),RXYCL) ,RYX(L) 
WRibeeee, 200) RCI) .RY(L),RXYCL),RYXCL) 
WRITE (9,201) 

CONTINUE 

WRITE (9,201) 

FORMAT (2X,4(2X,F14. 9)) 

FORMAT (' '‘) 


sesevevesesesesevesesevesevesesesese cece sesevesesesevesevevesevedslesesesesessseskckseslesvcsesleskeskskslesleskotsescoksesesksesssroks'esk 
INPUT THE DESIRED AR AND MA ORDERS THEN DEFINE INITIAL CONDITIONS 


WRITE (6,3) 

FORMAT (/' ENTER THE DESIRED AR ORDER: ') 
READ (6,*) NS 

WRITE (6,4) 

FORMAT (/' ENTER THE DESIRED MA ORDER: ') 
READ (6,*) NT 


KS=0 

KT=0 
A(0)=1. 0 
VX=RX(0) 
VY=RY(0) 
VXY=-RYX(0) 
VG=RX(0) 
VZ=RY(0) 
VGZ=-RYX(0) 
TXB(0)=1. 0 
TyACOj=1) 0 
GBCO)=1.6 
ZTA(0)=1. 0 


3, 


(C2 Ga 1G) Ca 


301 


302 


303 


ae Ee aD EG i) 
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Sesescacdedevevovesevevcvesevevevevesesevese sevevevesese se dese sevese sieve sevetove te Jevese seskeseve Tess vese se se se tese sete se 1086 te ee ae Ie 
ESTIMATE THE ARMA PARGMETERS 


IF CNIVEOQ? OVAND No ate NS) teen 
KS=KS+1 
CALL NEWARCKS,KT,VX,VY,VXY,VG,VZ,VGZ,TXA, TXB, TYA, TYB,GA,GB,ZTA,ZTB 
GRA SRY kay  RYAG A,B) 
GOTO 300 
ELSE 
TP (NS. EQ. 0. AND. Kia NT) Pa 
KT=KT+1 
CALL NEWMACKS ,KT,VX,VY,VXY,VG,VZ,VGZ,TXA,TXB,TYA,TYB,GA,GB,ZTA, ZTB 
GRA GRY. RAL, Rin ASD) 
COTO 201 
ENDIF 
ENDIP 


IF (NS.NE.0.OR.NT.NE. 0) THEN 
IF (NS. GE. NT. AND. NT. NE. 0. AND. KT. LT. NT) THEN 
KS=KS+1 
CALL NEWAR(KS,KT,VX,VY,VXY,VG,VZ,VGZ,TXA, TXB, TYA, TYB,GA,GB, ZTA, ZTB 
& ,RX,RY,RXY,RYX,A,B) 
KT=KT+1 
CALL NEWMACKS,KT,VX,VY,VXY,VG,VZ,VGZ,TXA,TXB, TYA, 1YB,GA,GB, 2) aayame 
&,RX,RY,RXY,RYX,A,B) 
GOTO 302 
ELSE 
IF (KS.LT.NS) THEN 
KS=KS+1 
CALL NEWAR(KS ,KT,VX,VY,VXY,VG,VZ,VGZ,TXA, TXB, TYA, TYB,GA,GB, ZTA, ZTB 
& ,RX,RY,RXY,RYX,A,B) 
GOTO 303 
ENDIF 
ENDIF 
ENDIF 


sevesevedolcdedsdedcdelcdedesvedevcdedevcledsdedcdededcdedevededcvevcdsvedevedeveds ievededevevevs Tevet veve ke Tevet de voce de He v6 Hee 
PRINT ESTIMATED ARMA PARAMETERS 


WRITE (*,211) 

WRITE (9,211) 

WRITE, (75210) (AGK)) Kade) 
WRITE €9.210)) (ACK)... Kelp es 
WRITE (*,211) 

WRITE (9,211) 

FORMAT(' ') 

WRITE (*.210)°@SCK) K-01 


WRITE (GomeOn) (BCK). K=0,KT) 
Ze HORMAT (Game 1 a( 2X 13. 10)) 

STOP 

END 


Sevesvoedcscdeseseveslcscdkevevevededevedesedevesescvededesesesdevedesesevevedesk dledeslevesesicvesesededesedesevevesesesesedesvvesodevevs 


SUBROUTINE TO COMPUTE CORRELATION TERMS 


ee al ae Ma > Hs Se os I 


SUBROUTINE CORREL(N,LAG,X,Y,RX,RY,RXY,RYX,NK1) 
REAL X(0:NK1),Y(0: NK1) ,RX(0: 2000) ,RY(0: 2000) ,RXY(0: 2000) 
REAL RYX(0: 2000) ,SUM1,SUM2 ,SUM3,SUM4 
DO 70 K=0,LAG 
NJ=N-1-K 
SUM1=0 
SUM2=0 
SUM3=0 
SUM4=0 
ANK=NJ 
DO 60 J=0,NJ 
SUM1=SUM1+X(J+K)*X(J) 
SUM2=SUM2+Y( J+K)*Y(J) 
SUM3=SUM3+X(J+K)*Y(J) 
SUM4=SUM4+X(J)*¥(J+K) 
60 CONTINUE 
RX(K)=SUM1/ANK 
RY(K)=SUM2/ANK 
RYX(K)=SUM3/ANK 
RXY(K)=SUM4/ANK 
70 CONTINUE 
RETURN 
END 
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APPENDIX B. SUBROUTINE FOR MAIN PROGRAM 


SUBROUTINE NEWAR(KS,KT,VX,VY, VXY,VG,VZ,VGZ,TXA, TXB, TYA,TYB,GA,GB,Z 
&TA,ZTB,RX,RY,RXY,RYX,A,B) 


THIS SUBROUTINE COMPUTES AR PARAMETER VALUES FOR AN ARMA MODEL AS 


THE AK ORDER INCKEASE Ss 5 YaGrr 


VARIABLE DEFINITIONS 


KS - CURRENT AR ORDER 

KT - CURRENT MA ORDER 

RX - AUTOCORRELATION DATA OF INPUT VN 

RY ~ AUTOCORRELATION DATA OF OUTPUT Y 

RXY - CROSSCORRELATION DATA OF INPUT AND OUTPUT 

RAY - CROSSCORRELATION DATA OF OUTPUT AND eiNEU 

TXA - VECTOR CONTAINING AR COEFFICIENTS FOR FORWARD PREDICTION 
OF INPud: 

TXB - VECTGR CONTAINING MA COEFFICIENTS FOR FORWARD PREDICTION 
OF INPUT 

TYA - VECTOR CONTAINING AR COEFFICIENTS FOR FORWARD PREDICTION 
OF OUTPUT 

TYB - VECTOR CONTAINING MA COEFFICIENTS FOR FORWARD PREDICTION 
OF QUTEUT 

GA - VECTOR CONTAINING AR COEFFICIENTS FOR BACKWARD PREDICTION 
OF INEUD 

GB - VECTOR CONTAINING MA COEFFICIENTS FOR BACKWARD PREDICTION 
Ob INPUl 

ZTA - VECTOR CONTAINING AR COEFFICIENTS FOR BACKWARD PREDICTION 
OF OUTPUT 

ZTB - VECTOR CONTAINING MA COEFFICIENTS FOR BACKWARD PREDICTION 
OF OU Pum 

C ARRAY WHICH STORES CURRENT VALUES OF TXA 

D ARRAY WHICH STORES CURRENT VALUES OF TYA 

E ARRAY WHICH STORES CURRENT VALUES OF GA 

iy ARRAY WHICH STORES CURRENT VALUES OF ZTA 

Is ARRAY WHICH STORES CURRENT VALUES OF TXB 

Q ARRAY WHICH STORES CURRENT VALUES OF TYB 

R ARRAY WHICH STORES CURRENT VALUES OF GB 

S) ARRAY WHICH STORES CURRENT VALUES OF ZTB 

A - VECTOR UPDATED AR COEFFICIENTS OF ARMA MODEL 

B - VECTOR CONTAINING UPDATED MA COEFFICIENTS OF ARMA MODEL. 

TAU1 - CONSTANT COMPUTED FROM CORRELATION DATA AND PREDICTION 
ERROR COEFFICIENTS. 

TAU2 - CONSTANT COMPUTED FROM CORRELATION DATA AND PREDICTION 
ERROR COEFFICIENTS. 

TAU3 - CONSTANT COMPUTED FROM CORRELATION DATA AND PREDICTION 
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aro 
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TAU4 


XMU1 
YMU1 
MU2 
XMU3 
YMU3 
MU4 
MUS 
DET 


ERR 


ERROR GOEFEICIENTS. 
CONSTANT COMPUTED FROM CORRELATION DATA AND PREDICTION 
DeeOkecenel TOLENTS. 


CORE EICIEN® 
COBEEICLIENT 
COnmn ie Lind 
COMER ICIENGT 
COERELCrEN 
CORPFICIENG 
COEEEICIEN® 
DETERMINANT 
NOV Vay 


OF 
(OF 
OF 
Or 
OF 
Or 
OF 
OF 


AR-TYPE RECURSIVE FORMULA 
AR-TYPE RECURSIVE FORMULA 
AR-TYPE RECURSIVE FORMULA 
AR-TYPE RECURSIVE FORMULA 
AR=TYEE SRECURSIVE FORMULA 
AR-TYPE RECURSIVE FORMULA 
AR-TYPE RECURSIVE FORMULA 
PREDICTION ERROR MATRIX COMPOSED OF 


ERROR BETWEEN REFERENCE MODEL OUTPUT AND LATTICE 
REALIZATION OUTPUT. 


DIMEN SONS ACO. 20000 kYC 0; Z000) ,_RXYCO: 2000) ,RYACO: 2000) ,AC 0: 21) 
aN See NCO 2h) PxACO; 21) TYACO: 21), 1TXB(0: 21), TYBC 0: 21) ,GAC0: 21) 
emer 0:2  ZmACO> 21. Z0BC Or 21) ,CC0: 21) (DCO: 21) ,EC 0: 21) 

Di ee Nee On 21) yOC0721),RC0: 21),50C0: 21), FCO: 21) 

REAL MU5 ,MU2 ,MU4 


COMPUTE VALUES FOR TAU1 THROUGH TAU4 


T1S=0 
T2S=0 
T3S=0 
T4S=0 


KI=KS-1 


DO 10 I=0,KI 
T1S=T1S-RYX( I+1)*ZTA(KI-I) 
T2S=T2S+RY(I+1)*ZTA(KI-I) 
T3S=T3S+RY( I+1)*GA( I) 
T4S=T4S+RY(I+1)*ZTA(1) 

CONTINUE 


T1T=0 
T2T=0 
T3T=0 
T4T=0 


DO 20 J=0,KT 
T1IT=T1T+RX( J+1)*ZTBCKT-J) 
P2t 02 ti CI+ 1) ZT BK) 
Pol tetek i XC Ki -Kitits) 6nd) 
T4T=T4T-RYX( KI -KT+1+J)*ZTB( J) 
CONTINUE 
TAU1=T1S+T1T 
TAU2=T2S+T2T 
TAU3=T35S+T3T 
TAUS=T4S+T4T 
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COMPUTE VALUES FOR THE REFLECTION COEFFICIENTS 


XMU1=-TAU1/VZ 

YMUI=-TAUZ/ VZ 

MU2 =-VG2Z/VZ 

DET =VX*VY-VXY*VXY 

XMU3=-( VY*TAU1-VXY*TAU2) /DET 

YMU3S=~-( ~VAY*TAUIFV 4 FA02 7 DED 

MU& =(VGZ*TAU4S-VZ*TAU3) /C( VGZ*VGZ-VG*VZ) 
MUS =MU4*MU2 


COMPUTE ARMA COEFFICIENTS 


DO 16 K=0,KS 

C(K)=TXA(K) 

D(K)=TYA(K) 

E(K)=GA(K) 

F(K)=ZTA(K) 

CONTINUE 

DO 45 J=1,KS 

TXA( J)=C( J)+XMU1*F(KS-J) 
TYA(J)=D(J)+YMU1*F(KS-J) 
GA( J)=E( J-1)+MU2*F(J-1) 

ZTA( J)=F( J) +XMU3*C( KS-J) +YMU3*D( KS-J)+MU4*E( J-1)+MUS*F(J-1) 
CONTINUE 

DO 31 K=0,KT 

P(K)=TXB(K) 

QCK)=TYBCK) 

S(K)=ZTB(K) 

R(K)=GB(K) 

CONTINUE 

DO 55 J=1,KT 
TXB(J)=P(J)+XMU1*S(KT+1-J) 
TYB( J)=Q( J)+YMU1*S(KT+1-J) 
GB(J) =R(J)+MU2*S(J) 
ZTB( J)=S( J+1)+XMU3*"*P( KT-J)+YMU3*Q( KT-J)+MU4*R( J) +MU5S*S(J) 
CONTINUE 

WRITE (*,176) KS 

WRITE (9,176) KS 

FORMAT( 12) 

WRITE (*,175) (ZTB(K), K=0,KT) 

WRITE (9,175) (ZTB(K), K=0,KT) 
FORMAT (4(1X,F10.5)) 


UPDATE ERRORS 
FORMAT (/' S UPDATE ERROR IS: ',F15. 10) 
VX =VX+XMU1*TAU1 


VY =VY+YMU1*TAU2 
VXY=VXY+XMU1*TAU2 
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70 


VG =VG+MU2*VGZ 

VZ=VZ+( XMU3*TAU1+YMU3*TAU2 )+MU4*TAU3+MU5*TAUG 
VGZ=TAU3+MU2*TAUS 

ERR=VY -( VXY"2) /VX 

WRITE (**,650) ERR 

WRITE (9,650) ERR 

WRITE (*,888) VX,VY,VG,VZ 

FORMAT (4(1X,F10. 6)) 


COMPUTE MODEL COEFFICIENTS 


DO 65 J=1,KS 
A( J)=TYA( J) -TXA( J) *VXY/VX 
CONTINUE 
DO 70 J=1,KT 
B( J)=TYB( J) -TXB( J)*VXY/VX 
CONTINUE 
B(0)=-VXY/VX 
RETURN 
END 


ag 
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APPENDIX C. SUBROUTINE FOR MAIN PROGRAM 


SUBROUTINE NEWMA(KS,KT,VX,VY,VXY,VG,VZ,VGZ,TXA,TXB,TYA,TYB,GA,GB,Z 
&TA,ZTB,RX,RY,RXY,RYX,A,B) 


THIS SUBROUTINE COMPUTES MA PARAMETER VALUES FOR AN ARMA MODEL AS 
THE MA ORDER INCREASES BY ONE. 


VARIABLE DEFINITIONS 


KS - CURRENT AR ORDER 

KY - CURRENT MA ORDER 

RX - AUTOCORRELATION DATA OF INPUT VN 

RY - AUTOCORRELATION DATA OF OUTPUT Y 

RXY - CROSSCORRELATION DATA OF INPUT AND OUTPUT 

RXY - CROSSCORRELATION DATA OF OUTPUT AND INPUT 

TXA - VECTOR CONTAINING AR COEFFICIENTS FOR FORWARD PREDICTION 
OF ANEUE. 

cXB - VECTOR CONTAINING MA COEFFICIENTS FOR FORWARD PREDICTION 
OF ITNEUL 

TYA - VECTOR CONTAINING AR COEFFICIENTS FOR FORWARD PREDICTION 
OF OUTPUT 

(ier, - VECTOR CONTAINING MA COEFFICIENTS FOR FORWARD PREDICTION 
OF OULU 

GA - VECTOR CONTAINING AR COEFFICIENTS FOR BACKWARD PREDICTION 
OF INEUE 

GB - VECTOR CONTAINING MA COEFFICIENTS FOR BACKWARD PREDICTION 
OF ANEUT 

ZTA - VECTOR CONTAINING AR COEFFICIENTS FOR BACKWARD PREDICTION 
OF OUTPUT 

ZTB - VECTOR CONTAINING MA COEFFICIENTS FOR BACKWARD PREDICTION 
OF OUTEOL 

C - ARRAY WHICH STORES CURRENT VALUES OF TXA 

D - ARRAY WHICH STORES CURRENT VALUES OF TYA 

E - ARRAY WHICH STORES CURRENT VALUES OF GA 

7 - ARRAY WHICH STORES CURRENT VALUES OF ZTA 

iS - ARRAY WHICH STORES CURRENT VALUES OF TXB 

Q - ARRAY WHICH STORES CURRENT VALUES OF TYB 

R - ARRAY WHICH STORES CURRENT VALUES OF GB 

5 - ARRAY WHICH STORES CURRENT VALUES OF ZTB 

A - VECTOR UPDATED AR COEFFICIENTS OF ARMA MODEL 

B - VECTOR CONTAINING UPDATED MA COEFFICIENTS OF ARMA MODEL. 

TAUIP - CONSTANT COMPUTED FROM CORRELATION DATA AND PREDICTION 


ERROR GOEFETCTENTS, 

TAU2P - CONSTANT COMPUTED FROM CORRELATION DATA AND PREDICTION 
ERROR COEFFICIENTS. 

TAU3P - CONSTANT COMPUTED FROM CORRELATION DATA AND PREDICTION 
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ERROR COEFFICIENTS. 

TAUSP - CONSTANT COMPUTED FROM CORRELATION DATA AND PREDICTION 
ERROR COEFFICIENTS. 

Reeaee = COMPETCTENT OP MA-TYPE RECURSIVE FORMULA 

YETA1 = COEFFICIENT OF MA-TYPE RECURSIVE FORMULA 

ETA2 - COEFFICIENT OF MA-TYPE RECURSIVE FORMULA 

XETA3 - COEFFICIENT OF MA-TYPE RECURSIVE FORMULA 

Pol oe CObrhrGiEN? OF MA-TYPE RECURSIVE FORMULA 

ETAS4 - COEFFICIENT OF MA-TYPE RECURSIVE FORMULA 

Poe cOnrr I CLENT OF MA-TYPE RECURSIVE FORMULA 


DET - DETERMINANT OF PREDICTION ERROR MATRIX COMPOSED OF 
Vay VAY. 
ERR - ERROR BETWEEN REFERENCE MODEL OUTPUT AND LATTICE 


DIMENSION RX(0: 2000) ,RY(0: 2000) , RXY( 0: 2000) ,RYX(0: 2000) ,A(0: 21) 
DIMENSION B(0: 21),TXA(0: 21) ,TXB(O: 21) ,TYA(O: 21) , TYB( 0: 21) ,GA(0: 21) 
DIMENSION GB(0: 21), ZTA(0: 21) ,ZTB(0: 21) ,C(0: 21) ,D(0: 21) 

DIMENSION E(0: 21),F(0: 21),P(0: 21) ,Q(0: 21) ,R(0: 21) ,S(0: 21) 


COMPUTE VALUES FOR TAU1 PRIME THROUGH TAU4S PRIME 


T1TP=0 

T2EP-0 

T3TP=0 

T4TP=0 

KJ=kT-1 

DO 10 I=0,KJ 
TITP—TITPeRXCi+1)*6B(KI=1) 
a2 T= 7 sGBCKI=1 ) 
ete reark Ai) ZU BCL) 
T4TP=T4TP+RX( I+1)*GBC I) 
CONTINUE 

T1iSP=0 

T2SP=0 

T3SP=0 

T45P=0 

DOR ZO I—02KS 
GiSe=TtSr shy <cd+))*+GA(KS-J) 
T2SP=T2SP+RY( J+1)*GACKS-J) 
tooP=(osPekny( Ki-kSt1+))*ZTAl J) 
T4SP=T4SP-RXYCKJ-KS+1+J)*GA( J) 
CONTINUE 

TAUIP=fTIP+TisSP 
TAU2P=T2TP+T2SP 
PAUSP=13EP+T3SP 
TAU4P=T4TP+T4SP 


COMPUTE VALUES FOR REFLECTION COEFFICIENTS 


XETA1=-TAU1P/VG 


61 


> 


16 
165 


41 
I7> 


LS el ee Ta Se aS 


66 


YETA1=-TAU2P/VG 

ETA2 =-VGZ/VG 

DET =VX*VY-VXY*VXY 

XETA3=-( VY *DAUIP-VAY*TAU2Z Py PED 

YETA3=-( -VAY*TAUIP+VX*TAU2P ) /DET 

ETA4S4 =(VGZ*TAUSP-VG*TAU3P ) /( VGZ*VGZ-VG*VZ) 
ETAS =ETA4S*ETA2 


COMPUTE ARMA PARAMETERS 


DO 16 K=0,KT 

P(K)=TXB(K) 

QCK)=TYB(K) 

R(K)=GB(K) 

S(K)=ZTB(K) 

CONTINUE 

FORMAT (4(1X,F10.5)) 

DO 45 J=1,kT 
TXB(J)=P(J)+XETA1*R(KT-J) 
TYB( J)=Q( J)+YETA1*R(KT-J) 
GB( J)=R( J) +XETA3*P( KT-J)+YETA3*Q(KT-J)+ETA4*S(J-1)+ETAS*R(J-1) 
ZTB( J)=S( J-1)+ETA2*R(J-1) 

CONTINUE 

DO 41 K=0,KS 

C(K)=TXA(K) 

D(K)=TYA(K) 

E(K)=GA(K) 

F(K)=ZTA(K) 

CONTINUE 

FORMAT (5(1X,F10.5)) 

DO 55 J=1,KS 
TXA( J)=C( J)+XETA1*E(KS+1-J) 
TYA( J)=D( J) +YETA1*E( KS+1-J) 
GA( J)=E( J+1)+XETA3*C( KS-J)+YETA3*D( KS-J)+ETA4*F(J)+ETAS*E(J) 
ZTA( J)=F( J) +ETA2**E( J) 

CONTINUE 


UPDATE ERRORS 


VX=VX+XKETA1*TAUI1P 

VY=VY+YETA1*TAU2P 

VXY=VXY+XETA1*TAU2P 

VG=VG+( TAU1P*XETA3+TAU2P* YETA3) +E TAG*TAU3P+ETA5S*TAUGP 
VZ=VZ+ETA2*VGZ 

VGZ=TAU3P+ETA2*TAUGP 

ERR=VY -( VXY****2) /VX 

WRITE (*,66) ERR 

WRITE (9,66) ERR 

FORMAT (/' T UPDATE ERROR IS: ‘',F15.10) 
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WRITE (*,889) VX,VY,VG,VZ 
FORMAT (4(1X,F10. 6)) 


COMPUTE MODEL COEFFICIENTS 


POmes J=1,KT 

BC J)=TYB( J) -TXBCJ)*VXY/VX 
CONTINUE 
DO 70 J=1,KS 

AC J)=TYAC J) -TXAC J) *VXY/VX 
CONTINUE 
BC O)=-VXY/VX 
RETURN 
END 
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APPENDIX D. ADAPTIVE LATTICE ALGORITHM PROGRAMNI 
THIS PROGRAM CALCULATES VALUES OF THE LATTICE COEFFICENTSAND 
OUTPUT OF AN ARMA DIGITAL LATTICE FILTER USING AN ADAPTIVE 
LATTICE ALGORITHM. 


VARIABLE DEFINITIONS 


X - ARRAY OF INPUT DATA, COMPUTER GENERATED WHITE GAUSSIAN 
NOISE WITH UNIT VARIANCE. 

AK = ARRAY OF DLAITICESCGEr relents. 

R * ARRAY OF LATTICE COERPICLIENTS: 

BO - TERMINAL CONDITION OF LATTICE REALIZATION. 

M - NUMBER OF LATTICE STAGES (EQUIVALENT TO ORDER OF ARMA 
rODEL). 

EXF - ARRAY OF FORWARD PREDICTION ERRORS FOR INPUT X. 

EXB - ARRAY OF BACKWARD PREDICTION ERRORS FOR INPUT X. 

EXBD - ARRAY OF DELAYED BACKWARD PREDICTION ERRORS FOR INPUT X. 

Bar - ARRAY OF FORWARD PREDICTION ERRORS FOR OUTPUT Y. 

EYB - ARRAY OF BACKWARD PREDICTION ERRORS FOR OUTPUT Y. 


EYBD - ARRAY OF DELAYED BACKWARD PREDICTION ERRORS FOR OUTPUT Y. 
ERROR - DIFFERENCE BETWEEN REFERENCE MODEL OUTPUT AND LATTICE 
REALIZATION OUTPUT. 


YE - ARRAYS CONTAINING LATTICE COEFFICIENT VALUES AT EACH 
ITERATION. 

MU - CONVERGENCE CONSTANT. 

RHO - WEIGHT GIVEN TO CUURENT POWER LEVEL AT EACH STAGE OF THE 


LATTICE STRUCTURE. 

SIGK - POWER LEVEL USED TO NORMALIZE CONVERGENCE CONSTANT WHEN 
UPDATING AK LATTICE COEFFICIENTS. 

SIGR - POWER LEVEL USED TO NORMALIZE CONVERGENCE CONSTANT WHEN 
UPDATING R LATTICE COEFFICIENTS. 

SIGB - POWER LEVEL USED TO NORMALIZE CONVERGENCE CONSTANT WHEN 
UPDATING TERMINAL CONDITION BO. 


DIMENSION EXF(10),EXB( 10) ,EXBD( 10) ,EYF( 10), EYB( 10) ,EYBD( 10) eee 
DIMENSION AK( 10), X(9900)}, ERROR( 9900}, V(12)0yYeuere a 

REAL MU 

SIGK=1. 

SIGR=1. 

SIGB=1. 


INITIALIZE ARRAYS 
DO 5 I=1,10 


EXF(1)=0 
EXB(1)=0 
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EXBD(1)=0 
EYF(1)=0 
EYB(1)=0 
EYBD(1I)=0 
R(1)=0 
AK(1)=0 
CONTINUE 

DO 6 I=1,9000 
X(1)=0 

ERROR( 1)=0 
CONTINUE 


ENTER VALUE OF THE CONVERGENCE CONSTANT MU AND VALUE OF RHO. 


M=2 
N=300 

WRITE (6,*) 'ENTER MU' 
READ(6,*) MU 

WRITE (6,*) 'ENTER RHO' 
READ (6,%*) RHO 


GENERATE WHITE GAUSSIAN NOISE 


ISIZE = N 
IX = 152255 
ISORT = 0 
MUL = 2 
DO 7 K= 1,ISIZE 
CALL SRND(IX,V,12,MUL, ISORT) 
XT=-6.0 
DO 8 I=1,12 
XT=XT+V(I) 
X(K)=XT 
CONTINUE 


COMPUTE OUTPUT OF REFERENCE MODEL FILTER AND LATTICE STRUCTURE 
TE NeGeNPUTE THE ERKOR. 


REFERENCE MODEL 


Y3=0 

Y2=0 

Y1=0 

X3=0 

X2=0 

X1=0 

BO=1. 

DO 100 I=1,N 
YF=X(1)-0. 8°X1+1. 78*X2+0. 89*Y1-0. 25*Y2 
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YF=0. 5*X(1)-0. 4*X1+. 89*X2+0. 89*Y1-0. 25*Y2 

YF=X(1)-2. 7*X1+3. 21*X2-1. 595%X341. 95*Y1-1. 62*Y2+0. 54*Y3 
YF=0. 5*X(1)-0. 95*X1+41. 33*X2-0. 979X341. 69*Y1-0. 962*Y2+0. 2*Y3 
YF=X(1)+0. 2*X1-0. 35*X241. 4*Y1-0. 85*Y2 

YF=0. 5*X(1)-0. 2*X14+0. 445*X2+1. OFY1-0. 94*Y2 


Y3=Y2 
Y2=Y1 
YI=YF 
X3=X2 
X2=X1 
X1=X( 1) 


LATIICE PURGE 


EXF(1)=X(1I) 

EXB(1)=X(I) 

DO 10 K=1,M 

EXF( K+1)=EXF(K)+R(K)*EXBD(K) -AK(K)*EYBD(K) 
EYF(M+1)=BO*EXF(M+1) 

DO 20 K=1,M 
EYF(M+1-K)=EYF(M+2-K)+R(M+1-K)*EXBD(M+1-K) -AK(M+1-K)*EYBD(M+1-K) 
EYB(1)=EYF(1) 

DO 30 K=1,M-1 

EXB( K+1)=EXBD( K)+R(K)*EXF(K) -R(K)*EYF(K) 
EYB( K+1)=EYBD(K)+AK( K)**EYF(K) -AK( K)*EXF(K) 
YL=EYF(1) 

ERROR(1)=YF-YL 

ERR=YF-YL 

CALL UPDATE (R,AK,EYBD,EXBD,ERR,MU,M,SIGK,SIGR, BO,RHO) 
CSB=EXF(M+1)*EXF(M+1) 

SIGB=RHO:*SIGB+( 1-RHO)*CSB 

BO=BO0+( MU/SIGB)**ERR*EXF(M+1) 

DO 40 K=1,M 

EXBD( K)=EXB(K) 

EYBD( K)=EYB(K) 

DO 50 J=1,M 

YE(J,1)=AK(J) 

YE( J+M,1)=R(J) 

FORMAT (2(1X,F10.6)) 

CONTINUE 


PRINT THE ERROR AND VALUES OF THE LATTICE COEFFICIENTS. 


WRITE (*,200) (ERROR(K), K=1,N,10) 
WRITE (9,200) (ERROR(K), K=1,N,10) 
FORMAT (5(1X,F10. 6)) 

WRITE (9,209) 

FORMAT(' ') 

WRITE (*,201) (R(K), K=1,M) 
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WRITE (*,201) (AK(K), K=1,M) 
WRITE (9,201) (R(K), K=1,M) 
WRITE (9,201) (AK(K), K=1,M) 
WRITE (*,205) BO 

WRITE (9,205) BO 

FORMAT (F10. 6) 

FORMAT (5(1X,F10.6)) 


CALL PLOTTING ROUTINES TO PLOT ERROR AND LATTICE COEFFICIENTS. 


CALL PLOT (ERROR,N) 
CALL PLOT1 (YE,N) 
STOP 

END 


SUBROUTINE WHICH UPDATES LATTICE COEFFICIENTS. 


SUBROUTINE UPDATE(R,AK,EYBD,EXBD,ERR,MU,M,SIGK,SIGR, BO, RHO) 
DIMENSION RC 10) ,AK(10) ,EYBD( 10) ,EXBD( 10) 
REAL MU 

CSK=0. 

CSR=0. 

DO 20 J=1,M 
CSK=CSK+EYBD( J) *EYBD( J) *BO**2 
CSR=CSR+EXBD( J)*EXBD( J) *BO**2 
SIGK=RHO*SIGK+( 1-RHO)*CSK 
SIGR=RHO*S IGR+( 1-RHO)*CSR 

DO 10 J=1,M 

RC J)=RCJ)+CMU/SIGR)*ERR*EXBD( J)*BO 

AKC J)=AK( J) -CMU/SIGK )*ERR*EYBD( J )*BO 
CONTINUE 

RETURN 

END 


Pood ice oUNE Oe LOT ERROR 


SUBROUTINE PLOT(Y,N) 
DIMENSION Y(N),X( 9900) 
DO 10 J=1,N 

Sd) 

CALL TEK618 

CALL PRTPLT(72,6) 
CALEGSUERPAC ADAPTIVE , A’ 53) 
CALL ORESHOC ALL) 

CALL PAGE(8. 50,6. 0) 
CALL HWROT( ' AUTO' ) 
CALI XIN DAX 
CALLEAREA2D( 5510. 3.0) 
CALL HEIGHT(0. 14) 

CALL COMPLX 
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CALL SHDCHR(90: 0,1,0.002) 

CALL HEADIN( 'LEARNING CURVESS' ,100,2.0,1) 
CALL XNAME( 'ITERATIONSS' ,100) 

CALL YNAME('ERRORS' ,100) 

CALL MESSAG(' ADAPTIVE FILTER  §$!',100,3.0,-0.8) 
CALL THKFRM(0. 03) 

CALL FRAME 

CALL GRAF(0, ‘SCALE {No-SeCOmmsOAne. ca 00> 
CALL THKCRV(0. 02) 

CAT) CURVEC X. YeNe0> 

CALL ENDPL(0) 

CALL DONEPL 

RETURN 

END 


SUBROUTINE TO PLOT LATTICE COEFFICIENTS 


SUBROUTINE PLOT1 (YE,N) 

DIMENSION YE(6,9900) ,X(9900) ,Y¥( 9900) , YD(9900) ,A(10) 
TRUE VALUES OF THE PARAMETERS 

A(1)=-0. 240719 

AC2)=0.125 

A(3)=-0. 195719 

A(4)=0. 8900 

A(5)=0. 89 

DO 10 J=1,N 

X(Jj=J 

CALL TEK618 

GALE -PRIELIC/ 2.6) 

CALL SHERPA( 'MENNECKE’ ,'A' ,3) 

PRINT SHERPA FILE: SHERPA XXYYZZXX SHGRAPH A 
CALL RESET( ALL ) 

CALL PAGECS.50,.4cn0) 

CALL HWROT('AUTO' ) 

CALL MINTAX 

GALL AREAZDGS.0..3-10) 

CALL HEIGHT(0. 14) 

CALL COMPLX 

CALL SHDCHR(90.0,1,0.002,1) 

CALL HEADIN(' PARAMETERSS' ,100,2.0,1) 
CALL XNAME('ITERATIONSS' , 100) 

CALL YNAME('MAGNITUDES' , 100) 

CALL MESSAG('ADAPTIVE ARMA LATTICES’ ,100,3.0,-0. 8) 
CALL THKFRM(0. 03) 

CALL FRAME 

CALL GRAF(0,'SCALE',N,-1.00, SCALE = i12@) 
GALE THRERV C0. C2) 

TO PLOT ESTIMATES 

DOS200K—14 
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DO 30 J=1,N 
Y(J)=YE(K,J) 

CALL CURVE(X,Y,N,0) 
CONTINUE 

TO PLOT TRUE PARAMETERS 
DO 40 K=1,4 

DO 50 J=1,N 
Y(J)=A(K) 

CALL DASH 

CALL CURVE(X,Y,N,0) 
CONTINUE 

CALL ENDPL(0) 

CALL DONEPL 

RETURN 

END 
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APPENDIX E. DERIVATION OF DIFFERENCE EQUATION FOR TWO 
STAGE ARMA LATTICE FILTER 


The lattice filter is described by the expressions for the forward and backward pre- 


diction errors, equations (3.6) and (3.7) respectively, 


> 


ey; (k) =x x(k) + wy x(k — l)— wy 5 (k — 1) 

ef (k) = ef (A) + 9 ep (A — 1) — wy ep (K-21) 

ep (k) = x(k — 1) + wy x(A) — 9g (A) (E-1) 
ep (k) = (A — 1) — wy x(k) + 043 3A) 

ef (k) = ef, (kK) + 3 ey (A — W) — 5 (kK -1) 


and the expression for the filter output, 
y(k) = e7 (4) Fi oe I)- 3 v(kK — 1) (E — 2) 
Substituting for e} (A) in equation (E-2) yields, 
x ie sees 2 yy: DEY ge a 
MK) = 67, (A) + ge, (A — 1) — wy gp (A — 1) + wg x(k — I — yA-—1) (E—3) 
Substituting for e;, (A — 1) and e} (A — 1) in equation (E-3) 


(kh) = ej, (A) + Wy Be (kK — 2) + ro x(k —1)- Wy y(k — 1) | 
— wy f sk — 2 — wy x(k — D+ wp s(k-1) | (E —4) 
me x(k —l)—w, {k— 1) 


Subsututing for ef (k) in (E-4) we obtain 


Vk) = ep (A) + 15 a (k-—1)—- We 5, (A—I)+ Wi x(k — 2) + Wy ro x(k — 1) 
th nie 1) us ie = Oe Wy rede "5 "3 Wk — 1) (4S) 
ee x(k a 1) = Wy V{k a 1) 


From (E-1), substitute for e; (A — 1) and e, (A — 1) in (E-5) to obtain 
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y(k) =e (k) + wy [ x(k — 2) + wy x(k - I wk - 1) J 
- Wwe [ (hk —2)- Wy x(k —1) + 13 v(k — 1) ] 
- wi x(k —2)+ wi Ww, x(k —1)- wi Wy Wk — 1) (20) 
i555 y(K — 2) + Wy Ww, x(k —1)- Ws Wy y(K — 1) 
+ wy x(k — 1) — w, x(k — 1) 


Now, from (E-1), substituting for e7 (k) in (E-6) vields 


y(k) = x(k) + 15 x(k -—1)- W, yk — 1) + ws x(k —2)+ WS i x(k — 1) 
_ ws Wy y(K — 1) we y(K — 2) + we Wy x(k —1)- we Wy y(k — 1) 


Be, xe yee ws 15 x(k — 1) ah. i wi y(k — 2) et 
+ w5 wy x(k — 1) — wh wy p(k — 1) toy x(k — 1) — 5 yk - 1) 
Grouping terms we get, 
y(k) = x(A) + (069 + wy ep Hoey wy ey wg Fwy w+ wg) x(k — 1) 
+ (wy + wy) x(A — 2) a 


_ (wy =O, 5 ae re We ae Ww + Wy ws + 3) (Kk =) 
2 2 4p. 
— (wy + w3) y(K — 2) 


From the gradient estimator and coefficient update equations we know that the follow- 


ing relationships among lattice coefficients are true 


an 


] ] 2 ] ] 2 p 
Wp SWzp WPS Wy Wy Shy WSS Wy aa 9) 


Using these equalities in equation (E-8), we obtain the final expression for the difference 
equation. 


(hk) = x(k) + 2(w, + rw, we + wy w5) x(k — 1) + 2 WS x(k — 2) 


E— 10) 
= 2(w, + i5 ws + WY We y¥(K — 1) -2 wy THUS 7) | | 
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